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1  Summary 


Signal  classification  or  automatic  modulation  classification  is  an  area  of  research  that  has  been 
studied  for  many  years,  originally  motivated  by  military  applications  and  in  current  years 
motivated  by  the  development  of  cognitive  radios.  Its  functions  may  include  the  surveillance  of 
signals  of  interest  and  providing  information  to  blind  demodulation  systems. 

The  problem  of  classifying  Code  Division  Multiple  Access  (CDMA)  signals  in  the  presence 
of  Additive  White  Gaussian  Noise  (AGWN)  is  explored  using  Decision  Theory.  Prior  state- 
of-the-art  has  been  limited  to  single  channel  digital  signals  such  as  MPSK  and  QAM,  with 
few  limited  attempts  to  develop  a  CDMA  classifiers.  Such  classifiers  make  use  of  the  cyclic- 
correlation  spectrum  for  single  user  and  feature-based  neural  network  approach  for  multiple 
user  CDMA.  Other  approaches  have  focused  on  blind  detection,  which  could  be  used  for 
classification  in  an  indirect  manner. 

The  discussion  is  focused  on  the  development  of  classifiers  using  the  average  likelihood 
function.  This  approach  will  ensure  that  the  development  is  optimal  in  the  sense  of  minimizing 
the  error  in  classification  when  compared  with  any  other  types  of  classification  techniques. 
However,  this  approach  has  a  challenging  problem:  it  requires  averaging  over  many  unknown 
parameters  and  can  become  an  intractable  problem. 

This  research  was  successful  in  reducing  some  of  the  complexity  of  this  problem.  Starting 
with  the  definition  of  the  probability  of  the  code  matrix  and  the  development  of  the  likelihood 
of  MPSK  signals,  it  was  possible  to  find  an  analytical  solution  for  CDMA  signals  with  a 
small  code  length.  Averaging  over  matrices  with  the  lowest  Total  Squared  Correlation  (TSC) 
allowed  simplifying  the  equations  for  higher  code  lengths.  The  resulting  algorithm  was  tested 
using  Receiver  Operating  Characteristic  Curves  and  Accuracy  versus  Signal-to-Noise  Ratio 
(SNR).  The  algorithm  that  classifies  CDMA  in  terms  of  code  length  and  number  of  active 
users  was  extended  to  different  complex  types  of  CDMA  under  the  assumptions  of  full-loaded, 
underloaded,  balanced  and  unbalanced  CDMA,  for  orthogonal  or  quasi-orthogonal  codes,  and 
chip-level  synchronization. 
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2  Introduction 

2.1  Modulation  Classification 

As  the  meaning  suggests,  modulation  classification  is  the  determination  of  the  modulation 
type  of  unknown  signals  for  the  purpose  of  identification,  as  it  is  often  the  case  in  military 
applications,  or  for  channelizing  the  signal  through  automating  signal  processing  algorithm, 
as  it  is  the  case  of  cognitive  radio  applications.  Modulation  classification  is  part  of  a  broader 
problem  known  as  blind  or  uncooperative  demodulation  the  goal  of  which  is  the  extraction  the 
information  contents  of  unknown  signals. 

An  important  consideration  in  solving  this  problem  is  the  selection  of  the  approach.  The 
literature  presents  three  common  types  of  classification  methods:  ad-hoc,  feature -based  and 
decision  theoretic  classification.  Ad-hoc  methods  are  of  problematic  because  they  are  based 
on  intuition,  offering  no  guarantees  in  the  performance  of  the  classifier  and  thus  becoming 
questionable  for  the  development  of  robust  telecommunication  systems.  The  feature-based 
approaches  attempt  to  classify  signals  based  on  the  extraction  of  ad-hoc  features  by  using  clus¬ 
tering  algorithm  such  as  neural  networks.  The  method  tries  to  infer  the  underlying  probabilities 
of  each  class  during  a  training  session  and  then  assigns  classes  by  measuring  some  distance 
metric  to  the  clusters  in  a  test  session.  The  use  of  these  methods  becomes  relevant  in  problems 
where  the  system’s  stochastic  model  is  either  incomplete  or  too  complex  to  be  described  in 
mathematical  terms.  Feature  based  methods  often  provide  an  acceptable  performance;  however, 
they  present  two  subjective  problems:  the  selection  of  features  and  the  selection  of  data  sets 
for  training.  These  weaknesses  also  leave  the  uncertainty  whether  optimal  classification  is 
achieved  because  these  methods  often  deal  with  a  performance  that  depends  on  the  training 
data  set  selected.  A  third  method  for  classification  is  model-based  decision  theoretic.  The 
method  makes  use  of  the  Bayes  Criterion  for  developing  mathematical  rules  that  guarantees 
optimal  performance  in  noise,  i.e.,  rules  that  guarantee  the  lowest  error  in  classification.  The 
method  is  suitable  in  problems  where  models  are  available  and  have  low  complexity.  Its  main 
disadvantage  is  the  development  of  rules  due  to  the  mathematical  complexity,  especially  when 
problems  deal  with  many  unknown  variables. 

The  research  presented  herein  is  dedicated  to  the  study  of  classical  decision  theoretic 
approaches  for  the  classification  of  CDMA  signals.  The  motivation  of  this  study  includes:  1.  the 
need  for  developing  applications  that  identify  and  exploit  unknown  signals,  2.  the  lack  of  study 
done  in  applying  average  likelihood  techniques  to  CDMA,  and  3.  the  benefit  of  implementing 
optimal  classification  method  that  provides  reliable  classification  in  a  noisy  channel. 

2.2  Research  Objectives 

The  research  objective  is  to  detect  or  classify  CDMA  signals  in  the  presence  of  noise.  The 
scope  of  this  research  includes:  the  classification  of  CDMA  under  various  scenarios  such  as 
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fully-loaded,  underloaded,  balanced  and  unbalanced  CDMA.  It  also  covers  the  classification  of 
four  types  of  CDMA  that  arise  from  a  combination  of  BPSK  and  QPSK  signals.  For  this  study, 
a  decision  theoretic  approach  was  chosen  because  of  the  lack  of  methods  for  the  classification  of 
CDMA,  this  method  would  serve  as  a  future  benchmark  for  comparing  alternative  classification 
methods. 

An  initial  assessment  of  the  problem  made  evident  that  classification  of  CDMA  is  extremely 
challenging.  The  first  major  problem  is  trying  to  detect  CDMA  signals  that  exhibit  noise-like 
characteristics.  Such  detection  would  need  to  rely  on  the  ability  to  identify  distinctive  statistical 
features  between  CDMA  and  noise  signals.  A  second  problem  deals  with  the  complexity  of  the 
CDMA  model.  In  CDMA,  the  large  number  of  unknowns  turns  the  classical  decision  theory 
into  a  difficult  task.  In  the  pursuit  of  this  goal,  the  research  will  consider  the  development  of 
MPSK  decision  theoretic  classifier  as  a  model  to  follow. 

The  prior  state-of-the-art  in  modulation  classification  lacks  of  modulation  classification 
techniques  for  CDMA.  A  survey  on  modulation  classification  methods  shows  that  prior  decision 
theoretic  approaches  have  been  applied  to  single  user  digital  signals.  The  few  approaches  related 
to  CDMA  are  based  on  the  concept  of  cyclostationary  features.  Two  of  these  approaches  are 
limited  to  single  user  CDMA.  A  third  method  is  not  a  classification  method  in  itself,  but  a  blind 
demodulation  method  that  can  be  used  for  parameter  estimation  in  a  Generalized  Likelihood 
Ratio  Test  (GLRT)  classifier.  The  algorithm  has  been  tested  up  to  code  lengths  of  32  with 
8  active  users.  The  idea  of  performing  brute  force  demodulation  prior  to  the  detection  will 
be  ruled  out  as  an  alternative  due  to  its  overwhelming  computational  costs  and  convergence 
problems. 

The  classification  of  CDMA  problem  started  off  the  well-known  MPSK  classifier  proposed 
in  [4].  The  intuition  suggest  that  a  CDMA  classifier  is  some  sort  of  convoluted  form  of  an  BPSK 
classifier  because  CDMA  is  based  on  BPSK  signals.  This  reasoning  proved  to  be  somewhat 
correct  because  all  CDMA  classification  rules  result  in  a  BPSK  rule  when  the  number  of  active 
users  is  one.  The  development  of  MPSK  classifier  also  served  as  a  starting  point  to  other  optimal 
signal  classifiers  such  as  QAM  and  FSK. 

One  of  the  most  obvious  challenges  in  the  developing  an  Average  Likelihood  Classifier  for 
CDMA  is  averaging  over  unknown  variables.  The  number  of  unknowns  grows  proportional  to  the 
dimensions  of  the  code  matrix  used  for  generating  CDMA;  however,  typical  code  matrices  such 
as  Hadamard  matrices  occurs  at  code  lengths  of  powers  of  2.  For  these  codes,  the  dimensions 
grows  also  exponentially,  which  clearly  presents  a  problem. 

A  key  finding  in  this  research  was  the  proposition  of  a  weight  for  averaging  CDMA  codes. 
This  weighting  function  is  referred  in  this  discussion  as  the  probability  of  the  code  matrix.  This 
weight  is  based  on  the  Total  Squared  Correlation  and  was  a  key  factor  in  the  development  of  a 
simplified  decision  rule.  Other  consideration  is  the  form  of  the  simplified  average  likelihood 
function.  It  was  found  that  this  function  can  be  expressed  as  a  product  of  hyperbolic  cosine 
functions  and  decaying  exponentials.  Many  of  these  decaying  exponentials  are  a  function  of  an 
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introduced  parameter  referred  as  the  precision  of  the  probability  of  the  code.  Increasing  the 
precision  of  our  classifier  eliminated  many  terms  of  the  average  likelihood.  This  simplification 
is  important  because  the  likelihood  function  can  easily  become  intractable  or  show  numerical 
problems. 

The  major  contribution  of  this  study  is  the  development  of  the  first  average  likelihood 
classifier  for  multiuser  CDMA.  The  hypothesis  under  test  is  the  code  length  and  the  number 
of  users,  which  in  some  aspect  is  analogous  to  the  classification  of  M-ary  PSK  with  M  as 
the  hypothesis  under  test.  The  simplified  average  likelihood  of  CDMA  was  developed  using 
a  standard  procedure,  so  the  likelihood  function  can  be  used  against  any  other  signal  such 
as  MPSK  or  QAM.  In  this  research,  CDMA  was  tested  against  BPSK  and  QPSK,  which  are 
special  cases  of  CDMA  under  the  assumption  the  the  code  length  equals  one.  The  same  average 
likelihood  for  CDMA  made  of  BPSK  symbols  can  be  reused  and  applied  to  all  types  of  CDMA 
generated  from  combinations  of  BPSK  and  QPSK  symbols  without  adding  more  complexity. 

2.3  Roadmap 

The  research  path  is  provided  in  Figure  1.  It  shows  the  main  step  taken  to  solve  the  problem. 
The  objective  is  to  detect/classify  CDMA  signals  in  AWGN,  Box  A  in  the  diagram.  The  study 
begins  with  a  brief  discussion  of  Detection  Theory  and  how  this  theory  was  applied  to  the 
previous  state-of-the-art  as  discussed  in  Chapter  3.  In  the  case  of  single  user  signals,  modulation 
classification  associates  a  digital  signal  to  a  constellation.  The  hypothesis  is  the  signal  type  and 
in  the  case  of  MPSK  and  QAM  signals,  the  type  is  associated  to  the  size  of  the  constellation.  In 
the  case  of  CDMA,  the  modulation  type  will  be  associated  to  the  code  length  and  the  number  of 
active  users. 

The  based  model  for  developing  and  simplifying  the  average  likelihood  function  of  CDMA 
signals  is  the  BPSK,  Box  B  in  the  diagram.  Both  signals  are  based  on  binary  antipodal  variables 
and,  at  the  end  of  the  development,  their  likelihoods  share  some  similarities.  The  likelihood 
for  one  BPSK  symbol  is  very  simple.  It  is  expressed  as  a  product  of  an  exponential  function 
in  terms  of  the  signal-to-noise  (SNR)  ratio.  However,  the  likelihood  of  a  CDMA  is  extremely 
complicated  for  a  CDMA.  This  problem  was  solved  by  using  symbolic  algebra  algorithms  (Box 
C)  for  small  code  lengths  as  discussed  in  Chapter  4. 

After  analysing  code  lengths  of  2,  3  and  4,  the  research  turns  its  attention  to  the  development 
of  a  mathematical  formulation  CDMA  (Box  D)  with  the  following  assumptions: 

■  chip-synchronous  CDMA, 

■  frame-asynchronous  CDMA, 

■  type  1  (BPSK-code/BPSK-data), 

■  unknown  orthogonal  code  matrix  (unknown  permutations), 
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Fig.  1  Research  Path 
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■  unknown  data  vector, 

■  unknown  code  length  and  number  of  users  (hypothesis), 

■  fully-loaded. 

The  case  of  unbalanced  CDMA  and  averaging  over  Gaussian  Distributed  amplitude  devia¬ 
tions  over  a  nominal  energy  per  symbol  (Box  F).  This  work  has  been  presented  in  MILCOM 
2014.  [5]  The  same  development  was  slightly  modified  and  applied  to  other  types  of  CDMA 
generated  from  a  combination  of  QPSK  and  BPSK  symbols  (Box  E).  A  key  approach  was  to 
represent  complex  models  in  terms  of  real  block  matrices  and  apply  an  approach  similar  to  type 
1  CDMA.  A  technical  paper  has  been  submitted  to  MILCOM  2016.  [6] 

2.4  Organization 

This  publication  is  organized  in  six  sections.  This  introduction  covers  the  problem  of  modulation 
classification  and  discusses  briefly  the  research  objectives.  In  section  3,  the  discussion  covers 
some  technical  background  with  the  prior  state-of-the-art;  discusses  some  essential  concepts 
found  in  decision  theory;  introduces  the  potential  approaches;  and  provides  and  insight  of 
how  decision  theory  is  applied  to  the  classification  of  MPSK  signals,  section  4,  considers  the 
classification  problem  applied  to  CDMA  and  provide  the  unique  mathematical  framework.  The 
average  likelihood  is  simplified  in  section  5.  In  section  6  the  concept  is  extended  to  other  types 
of  CDMA  derived  from  combinations  of  BPSK  and  QPSK  symbols.  Finally,  section  7  provides 
a  summary  of  the  findings,  future  work  and  conclusions. 
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3  Modulation  Classification  Methods 


This  section  provides  a  brief  review  of  the  concepts  found  in  decision  theory  and  how  they 
are  applied  to  the  problem  of  modulation  classification.  Existing  methods  in  modulation 
classification  can  be  grouped  into  three  main  categories  which  are  discussed  in  the  following 
sections.  Feature -based  and  decision  theoretic  provide  a  strong  mathematical  framework  and 
become  the  main  subjects  in  the  discussion.  The  decision  theoretic  method  will  be  the  method 
of  choice  for  developing  CDMA  classification  rules  after  evaluating  the  prior  state-of-the-art 
and  their  weaknesses.  Selecting  a  decision  theoretic  approach  is  risky  research  due  to  the 
complexities  in  averaging  the  CDMA  model,  but  this  risk  is  compensated  by  developing  a  new 
algorithm  based  on  optimization  principles. 

3.0.1  Ad  Hoc 

As  mentioned  before,  Ad  Hoc  classifiers  consist  of  deriving  intuitive  rules  for  classification. 
Several  of  these  methods  have  been  discussed  in  [7]  and  applied  to  military  applications. 
Some  of  them  are  simple  as:  thresholding  signals,  counting  samples,  generating  histograms 
[8]  or  applying  mathematical  rules  such  and  the  well-known  Mth -Power  Classifier  for  MPSK 
signals.  The  construction  of  these  rules  do  not  allow  for  analytically  predicting  the  performance. 
Although  simple  rules  such  as  the  Mth  power  law  could  be  approximations  of  optimal  rules, 
they  are  commonly  considered  of  limited  significance  unless  their  performance  is  supported  by 
some  theoretical  framework. 

3.0.2  Feature  Based  Methods 

A  more  scientific  approach  consists  of  identifying  important  features  of  a  signal  and  try  to 
cluster  the  features  in  a  multidimensional  space.  The  different  clusters  will  be  associated  to 
a  modulation  type  by  using  a  clustering  technique  such  as  neural  network.  The  approach  is 
a  good  choice  when  the  there  is  little  known  about  the  stochastic  model.  Examples  of  these 
methods  are  the  Statistical-Moment  Based  Classifier  presented  by  [9].  The  rationale  is  that 
statistical  moments  derived  from  signal  parameters  such  as  frequency,  amplitude  or  phase  can 
provide  useful  classification  features  for  a  neural  network  classifier. 

A  reference  of  this  technique  is  the  method  proposed  by  [1]  for  classifying  single  user 
CDMA  using  cyclostationary  features,  i.e.,  the  peaks  found  in  the  cyclic  correlation  spectrum. 
These  features  are  obtained  by  applying  a  Fast  Fourier  Transform  (FFT)  and  averaging  over 
frequencies  as  illustrated  in  Figure  2.  Frequency  smoothing  is  used  to  find  significant  peaks 
which  are  fed  into  a  neural  network. 

The  selection  of  peaks  as  features  is  motivated  by  the  signal’s  cyclostationary  properties 
of  a  single  CDMA  user.  The  referenced  method  is  limited  to  a  single  spreading  sequence  and 
classifying  multiuser  CDMA  signals  is  out  of  the  scope  of  this  particular  development.  Under 
the  proposed  example,  the  sequence  resulting  from  the  modulation  of  a  QPSK  spreading  and 
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a  BPSK  modulation  is  another  QPSK  signal.  So  in  reality  this  neural  network  classifier  is 
processing  periodic  QPSK  signals  for  ultrawide  band  radar  applications  rather  than  for  multi¬ 
user  CDMA  for  telecommunications.  The  report  tested  the  classification  of  two  pseudo  noise 
sequences  with  code  lengths  15  and  31.  The  classifier  was  able  to  achieve  accuracies  between 
86.5  and  100  percent  for  SNR  —  —MB  and  SNR  =  °°  respectively. 


BPSK 


{L,U} 

{15,1} 

{31,1} 


L  =  code  length 
U  =  active  users 

Fig.  2  Feature-Based  Neural  Network  Classifier 

The  performance  of  a  neural  network  depends  on  these  features  and  so  there  is  no  guarantee 
that  the  features  selected  is  a  fairly  complete  representation  of  all  the  features  needed  for  a 
correct  classification.  In  addition  to  this  fact,  the  performance  of  clustering  methods  is  known 
to  depend  on  the  training  data. 

3.0.3  Decision  Theoretic  Methods 

Decision  theoretic  methods  deal  with  the  problem  of  optimization  of  metrics:  the  minimization 
of  a  cost  function  or  the  maximization  of  the  accuracy.  The  theory  is  based  on  classical  detection 
theory  [10].  Bayes  Theorem  provides  the  statistical  model  that  characterizes  the  source,  the 
channel  and  the  observation  space  as  shown  in  Figure  4.  In  our  problem,  the  source  generates  a 
finite  set  of  known  classes  { representing  the  modulation  types.  Each  class  is  associated  to 
probability  referred  as  the  prior  probability.  The  source  generates  an  element  of  a  class 

and  sends  it  through  a  noisy  channel  which  is  characterized  by  a  likelihood  probability 
or  the  probability  of  the  observation  when  a  class  in  a  noisy  channel.  The  observation  r  is  fed 
into  a  detector  that  contains  a  classification  rule  based  on  optimality  principles. 

Table  1  shows  four  types  of  outcomes  may  occur  in  a  binary  classification.  Similar  concepts 
can  be  extended  to  multiple  classes  by  constructing  an  error  probability  matrix  of  multiple 
dimensions.  The  sum  of  non-diagonal  elements  represents  the  total  probability  of  error  in 
classification  while  the  sum  of  diagonal  elements  represents  the  accuracy.  The  construction  of 
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(X0,XU... 

Source 

Class 

Channel 

Observation  ^ 

Decision 

} 

pCK) 

- ► 

K 

p(R\M) 

R 

Rule 

Fig.  3  Bayesian  Model 
Table  1  Error  Probability  Matrix 


True  Class: 

Class 

Class  M\ 

Decide  for 

True  Negative 
P(Decide 

False  Negative 
P(Decide 

Decide  for 

False  Positive 
P(Decide  M\  J^o) 

True  Positive 
P(Decide  \d%\) 

an  optimal  decision  rule  shall  maximize  the  sum  of  the  diagonal  elements  and  minimize  the 
non-diagonal  elements. 

The  decision  process  requires  establishing  thresholds  or  decision  regions  in  the  observation 
space  of  r.  A  cost  function  can  be  assigned  to  each  decision,  but  these  costs  are  commonly 
ignored  in  classification  problems.  The  optimization  of  the  decision  regions  for  achieving 
maximum  accuracy  is  equivalent  to  maximizing  the  posterior  probability  according  to  (1). 

Decide  for  —  arg max  p(MJ\r) p(r) 

=  argmax  p{r\M’)p(M>) 

czu 

When  the  priors  are  not  available,  a  decision  is  taken  by  choosing  the  maximum  likelihood 
as  shown  in  (2).  Ignoring  the  priors  is  equivalent  to  assuming  that  all  prior  probabilities  are 
equal. 

Decide  for  =  argmax  p(r\M>)  (2) 

r~  f  o/P  OS/?  \ 

In  binary  classification,  the  same  criterion  can  be  expressed  in  a  form  of  a  Likelihood  Ratio 
Test  (LRT)  as  shown  in  (3).  The  determination  of  the  likelihood  function  and  its  simplifica¬ 
tion  is  an  important  step  in  deriving  optimal  rules  for  modulation  classification.  Computing 
the  likelihood  becomes  more  complicated  when  the  likelihood  involves  multiple  unknown 
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Table  2  Hypotheses  in  Various  Modulation  Classification  Problems 


Problem 

Hypothesis 

Cluster 

Detection  of  BPSK 

Jfe{±  1} 

point  in  the  C  plane 

Classification  of  MPSK 

M?e{M  =  2,4,...} 

constellation  in  the  C  plane 

Classification  of  CDMA 

jfe{L,u} 

vectors  a  G  ML 
aT a  ~L  U 

parameters. 

P(r\Jf\)  Decide^  p(j%) 

Pir  \J%)  Decide  Mi ]  P{M) 

Modulation  classification  deals  with  signals  with  unknown  parameters  in  noise.  The 
resulting  likelihood  is  conditioned  on  a  parameter  vector  0  that  is  unknown  to  the  classifier 
and  has  to  be  supplied  in  some  way  prior  to  classification.  Often,  the  parameters  are  treated  as 
random  variables  and  the  expectation  must  be  taken  over  all  the  possible  values  of  0  according 
to  (4)  in  order  to  remove  the  condition  according  to  the  Bayes  Theorem. 

p{r\Jif)  =  Eg{p(r|  JT,  0)}  (4) 

A  hypothesis  can  be  interpreted  as  the  label  of  a  region  belonging  to  a  class  in  the  observation 
space.  (See  Table  5)  In  detection  problems,  the  hypothesis  is  the  symbol.  In  blind  demodulation 
of  single  user  signals,  the  hypothesis  can  be  visualized  as  a  label  associated  to  a  constellation 
in  the  complex  plane.  In  the  case  of  CDMA,  the  hypothesis  can  be  visualized  as  a  label  to 
associated  to  a  collection  of  particular  vectors. 

3.1  Decision  Theory  for  Modulation  Classification 

Stochastic  processes  are  characterized  by  probabilities  that  evolve  over  parameters  such  as  time 
t  which  are  not  associated  to  a  random  variable.  A  stochastic  processes  n(t)  is  characterized 
by  the  statistical  moments.  Modulation  classification  often  deals  with  wide-sense  stationary 
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processes  characterized  by  having  a  zero  mean  and  constant  autocorrelation  function: 


E{n(t)}  =  0, 

Nn  (5) 

E{n(t)n*(t  +  At)}  =  ^8(At). 

Modulation  classification  also  deals  with  cyclostationary  processes,  i.e.,  processes  that  have  a 
periodic  nature  in  the  correlation  function  &  as  shown  in  (6).  The  parameter  a  in  this  equation 
represents  the  cyclic-correlation  period. 

3^  (At)  =  E{r(t)r*(t  +  At)} 

&(A  t  +  a)—  &(A  t)  ^ 


A  common  channel  found  in  telecommunications  is  the  Additive  White  Gaussian  channel  given 
by  (7).  The  modulated  signal  is  transmitted  through  a  noisy  channel  with  noise  n(t).  At  the 
receiver,  we  have  the  observation  in  the  form  of  a  process  r(t).  A  classification  rule  developed 
from  (1)  or  (4)  requires  using  random  variables.  It  would  be  impractical  to  characterize  the 
likelihood  using  an  infinite  number  of  observations  generated  by  r(t)  at  every  instant. 

r(t)  =  s(t,0)  +  n(t) 
n(t)  ~  <zK(jU,  a2) 


Parametrizing  the  likelihood  in  terms  of  t  is  inconvenient  because  the  observation  r(t)  has 
infinite  dimensions  and  the  likelihood  of  a  given  ensemble  of  r(t)  is  zero.  A  practical  way 
of  handling  processes  is  by  approximating  the  observation  to  a  finite  set  of  observations  by 
using  some  finite  approximation.  A  finite  set  of  observations  can  be  achieved  by  defining  a 
complete  orthonormal  (CON)  set  in  (8)  and  approximating  r(t)  in  terms  of  a  finite  set 

of  independent  coefficients  (r^(0)}. 


con  =  {y/k(t)}k=0:°° 

(Vi(t),Vj(t))  =  I  Vi(t)Vj(t)dt  =  8ij 

J  — oo 


In  modulated  sequences,  a  convenient  choice  for  our  basis  is  a  set  of  normalized,  non¬ 
overlapping  pulses  with  a  width  T  delayed  by  kT  according  to  (9). 


¥k(t) 


—  for  T  k  <  t  <  T(k  +  1)  and  k  Integer 
0  otherwise 


(9) 
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The  reduction  in  dimensionality  is  achieved  by  truncating  the  series  to  a  finite  number  of 
coefficients  ( K  <  °°)  according  to  (10).  The  likelihood  p ( { rk } K  M’)  will  be  constructed  from 
the  statistics  of  each  random  variables  rk.  For  a  linear  Gaussian  process  r(t),  the  statistics  of 
these  random  variables  are  also  Gaussian. 


n{0)  =  (r(t),}fk(t)) 
KO  ~  E  rk(0)Yk(t) 

k= 0 


(10) 


Eliminating  the  dependency  of  the  observation  parameters  requires  either  an  estimate  of  the 
parameter  vector  or  performing  the  expectation  over  all  the  unknowns  as  shown  in  (4).  If  the  last 
method  is  selected,  the  expectation  can  be  calculated  in  several  ways.  The  exact  computation  of 
the  expectation  results  in  average  likelihood  function  and  the  classifier  is  an  Average  Likelihood 
Ratio  Test  (ALRT).  An  approximation  of  the  expectation  results  in  Quasi- Average  Likelihood 
Ratio  Test  (QALRT).  The  substitution  of  0  by  estimated  parameters  results  in  a  Generalized 
Likelihood  Ratio  Test  (GLRT).  A  combination  of  any  of  these  methods  results  in  a  Hybrid 
Likelihood  Ratio  Test  (HLRT). 

3.2  Modulation  Classification  in  the  Literature 

Table  3  shows  a  list  of  known  methods  adapted  from  [2] .  Likelihood  functions  are  specific  to 
modulation  signals  and  their  respective  channel  models  which  may  be  seen  as  a  disadvantage 
because  changing  the  model  will  require  a  new  computation  of  the  likelihood  function.  The 
most  popular  classifiers  are  based  on  single  channel  BPSK,  QPSK,  MPSK  and  QAM  signals. 

Decision  theoretic  methods  are  less  abundant  in  the  literature  when  compared  to  feature- 
based  approaches.  This  can  be  attributed  to  the  complexity  in  the  development  of  such  algo¬ 
rithms.  Decision  theoretic  methods  based  on  likelihood  functions  require  good  signal  models 
and  few  unknown  parameters.  To  overcome  some  of  the  difficulties  in  calculating  the  average 
likelihood  function,  authors  recur  to  sub-optimal  approaches  such  as  QLRT,  GLRT  and  HLRT. 
The  case  of  CDMA  has  simple  models,  but  as  the  number  of  unknowns  increases,  the  averaging 
process  translates  into  higher  computational  costs. 

From  the  survey  of  decision  theoretical  methods,  it  is  evident  that  there  was  a  lack  of 
likelihood  methods  for  classifying  of  CDMA  signals  and  it  is  speculated  that  the  absence  of 
CDMA  methods  is  related  to  the  complexity  of  the  averaging  process.  A  contributions  of  this 
research  include  the  development  and  publication  of  classifier  for  BPSK-code/BPSK-signals  in 
2014  and  the  submission  of  a  similar  paper  on  complex  CDMA  classification  in  2016. 
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Table  3  Survey  of  Likelihood-Based  Classifiers 


Authors 

Classifier 

Modulation 

Unknowns 

Channel 

Sills 

ALRT 

BPSK,  QPSK, 
16QAM,  V29, 
32QAM, 

64QAM 

carrier  phase 

AWGN 

Wei,  Mendel 

ALRT 

16QAM,  V29 

- 

AWGN 

Kim,  Polydoros 

QLRT 

BPSK,  QPSK 

carrier  phase 

AWGN 

Sapiano,  Martin 

ALRT 

UW,  BPSK, 

QPSK,  8PSK, 
16PSK 

carrier  phase 

AWGN 

Long 

QLRT 

UW,  BPSK, 

QPSK,  8PSK 

carrier  phase, 
timing  offset 

AWGN 

Hong,  Ho 

ALRT 

BPSK,  QPSK 

symbol  level 

AWGN 

Beidas,  Weber 

QLRT 

32FSK,  64FSK 

phase  jitter 

AWGN 

Beidas,  Weber 

QLRT 

32FSK,  64FSK 

phase  jitter,  tim¬ 
ing 

AWGN 

Panagiotu 

GLRT,  HLRT 

16PSK, 

16QAM 

carrier  phase 

AWGN 

Chugg 

HLRT 

BPSK,  QPSK, 
OQPSK 

carrier  phase, 
signal  power, 
PSD 

AWGN 

Hong,  Ho 

HLRT 

BPSK,  QPSK 

signal  level 

AWGN 

Hong,  Ho 

HLRT 

BPSK,  QPSK 

angle  of  arrival 

AWGN 

Dobre 

HLRT 

BPSK,  QPSK, 
16QAM,  V29, 
32QAM, 

64QAM 

channel  ampli¬ 
tude,  phase 

flat  fading 

Abdi 

ALRT,  QLRT 

16QAM, 

32QAM, 

64QAM 

Channel  ampli¬ 
tude  and  phase 

flat  fading 

Vega-Irizarry, 

Fam 

ALRT 

CDMA  type  1 

code,  data  vec¬ 
tor 

AWGN 

Vega-Irizarry, 

Fam 

ALRT 

CDMA  types  2- 
4 

code,  data  vec¬ 
tor 

AWGN 

3.2.1  Generalized  Likelihood  Methods 

A  generalized  classification  method  for  CDMA  may  consist  of  making  a  estimate  of  the 
spreading  matrix  C  and  data  vector  b  with  code  length  L,  number  of  users  U  and  energy  per 
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symbol  E,  as  defined  in  the  CDMA  model  of  (11).  If  an  accurate  estimate  is  available,  a 
generalized  classification  rule  would  outperform  the  results  of  an  average  likelihood  classifier. 
However,  in  the  case  of  CDMA,  such  estimation  would  be  at  the  expense  of  high  computational 
costs  and  therefore  the  average  approach  would  seem  as  a  more  viable  approach.  For  comparison 
purposes,  we  consider  the  development  of  a  GLRT  using  Maximum  A  Posteriori  (MAP)  or 
Expectation  Maximization  for  estimating  unknown  parameters  prior  to  classification. 


Ce{±l}LxU 

9k  e  (±1}L 

nk~  jV  (0,Nq/2  I) 

for  k  =  0,  1 


(11) 


The  performance  of  generalized  methods  depends  on  the  accuracy  of  the  parameter  estimate, 
so  this  approach  can  be  wasteful  in  computational  resources  when  estimating  the  parameters  of 
wrong  classes.  For  example,  matrices  of  size  4x4  have  a  total  of  216  distinct  realizations  using 
BPSK  symbols.  An  algorithm  would  need  to  be  smart  enough  to  discard  many  useless  code 
matrices  that  would  never  be  used  in  CDMA  transmissions.  It  would  have  to  know  that  from 
a  total  of  216,  only  768  are  orthogonal  matrices  suitable  for  CDMA,  all  of  them  permutations 
of  one  single  matrix.  As  the  dimensions  of  the  code  matrix  increase,  this  kind  of  processing 
becomes  infeasible. 


3.2. 1.1  Maximum  Likelihood  Estimator 

The  parameters  C  and  b  can  be  estimated  using  a  maximum  likelihood  estimator.  Given  a 
likelihood  function  of  a  multivariate  Gaussian  stochastic  process  (12),  one  can  assume  the 
values  L  and  U  and  try  to  estimate  the  parameters.  Because  we  are  dealing  with  observations  jy. 
of  different  code  lengths,  this  would  require  reformatting  the  received  vector  in  different  code 
lengths. 


p({yt}K\^,C,  {b}K)  =  n  (aNl)L/2exP  (IP*  -  VEjicbt  if/MO  (12) 

Under  unknown  code  matrix  and  data  vector  samples,  the  problem  of  maximizing  the 
likelihood  of  a  CDMA  signal  (13)  does  not  have  a  closed  form  solution.  In  addition,  its  large 
number  of  variables  makes  the  computation  of  the  parameters  extremely  difficult  to  estimate  due 
to  the  limitation  of  numerical  methods  use  to  solve  the  equations.  Using  an  average  likelihood 
estimate  of  the  parameters  for  classification  purposes  will  not  solve  the  classification  problem 
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in  an  efficient  way. 


ViP(yk\LxUCDMA,C,b)  =  0 

^'£cH(yk-VE/LCbk)  =  0 
k 

Vcp(yk\LxUCDMA,C,bk)  =  0 
-^YM^-y/EjiCbk)  =  0 

k 


(13) 


3.2.1.2  Expectation  Maximization  Approach 

An  Expectation  Maximization  (EM)  approach  for  CDMA  for  blind  detection  [11]  provides  an 
iterative  way  of  estimating  the  code  matrix  C  as  an  updatable  parameter  and  the  data  vector  b  as 
a  hidden  random  variable  of  the  EM  algorithm.  EM  belongs  to  the  group  of  statistical  inference 
methods  in  which  the  parameter  is  adjusted  such  that  it  minimizes  the  distance  between  an 
empirical  distribution  and  a  model  distribution  [12]. 

The  proposed  algorithm  consists  of  four  major  steps. 

1.  Project  y  into  a  signal  space. 

2.  Define  El  as  the  updatable  code  matrix  given  by  (16). 

3.  Estimate  b  using  a  MMSE  estimator. 

4.  Maximize  the  expectation  to  obtain  an  update  of  £2. 

5.  Repeat  the  last  two  steps  until  convergence  is  achieved. 


The  algorithm  assumes  that  the  code  length  L  is  known  and  the  number  of  users  U  is 
determined  by  setting  a  threshold  on  the  singular  values  of  the  correlation  matrix  3%. 


&=[VS  V„ 


A,  0 
0  An 


[Vs  Vn]H  =  E{yf} 


(14) 


The  eigenvalues  of  the  signal  space  given  by  As  contain  the  energy  of  the  signal.  The 
eigenvalues  of  the  noise  space  given  by  A„  contain  only  noise  energy  as  shown  in  (15).  The 
eigenvectors  in  Vs  associated  with  the  signal  space  contain  information  about  the  code  matrix 
used  to  generate  the  CDMA  signal. 

Ei/L  +  No/2  for  signal  present 

A  r  •  ,  ’  \An)i,i 

0  for  noise  only 


0  for  signal  present 
No/I  for  noise  only 


(15) 
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The  observation  r  is  the  projection  of  y  into  the  signal  space.  This  is  an  estimate  value 
because  it  requires  a  threshold  to  separate  the  signal  from  the  noise  space  in  (14). 


r  =  VsHy 
Q  =  V”  C 


(16) 


The  updatable  matrix  Q  is  the  product  of  the  projection  and  the  code  C.  Both  Q  and  r  are 
parameters  of  the  likelihood  function  given  below. 


PpM l«>  = 


(17) 


The  data  vector  b  can  be  calculated  using  a  Minimum  Mean  Squared  Error  (MMSE)  detector 
[13]  as  shown  in  (18). 


b  =  sign{wT  r) 
w  —  arg  min  E  {b  —  b} 


(18) 


The  updated  Q  is  calculated  from  the  EM  algorithm  in  (19). 

a  =  argma  xp({rk}K,  {bk}K  |Q) 

a 


(19) 


The  computation  of  the  priors  p({bk}K)  presents  a  problem  in  the  algorithm  because 
computing  for  possible  2U  values  becomes  computationally  expensive.  Instead,  the  author 
has  decided  to  simplify  the  prior  probability  by  forcing  it  to  be  constant  p({bk}K)  —  1  /2UK. 

As  it  was  mentioned  before,  estimating  the  code  matrix  using  (19)  would  be  impractical  for 
modulation  classification  because  any  attempt  to  demodulate  a  wrong  hypothesis  results  in  a 
waste  of  computational  resources  and  processing  time.  Also,  the  blind  demodulation  algorithm 
has  been  tested  up  to  code  lengths  of  32  and  8  active  users.  Extending  the  algorithm  to  higher 
code  lengths  would  be  unusable  because  the  convergence  of  the  EM  would  degrade  as  the 
number  of  clusters  increases.  The  number  of  clusters  would  be  equivalent  to  the  number  of 
elements  in  the  code  matrix.  EM  is  a  method  of  choice  when  the  number  of  clusters  is  relatively 
small  for  quick  convergence,  which  is  not  the  case  for  classifying  higher  code  lengths.  So  in 
conclusion,  the  possibility  of  using  EM  blind  demodulation  method  for  constructing  GLRT 
modulation  classifier  can  be  ruled  out  as  a  viable  alternative. 
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3.2.2  The  Average  Likelihood  Method 

In  this  section,  the  average  likelihood  for  CDMA  approach  will  be  studied  by  following  the 
development  of  the  MPSK  classifier.  First,  it  would  be  assessed  if  constructing  a  simplified 
likelihood  function  is  feasible  and  meaningful  for  classification  purposes.  The  strategy,  if 
correctly  implemented  and  tested,  is  expected  to  reduce  to  a  BPSK  likelihood  function  for  code 
lengths  of  1 . 

3.2.2.1  MPSK  Model 


Ideal  MPSK  signals  in  (20)  are  constructed  from  a  sequence  of  orthonormal  pulses  (9)  with 
amplitude  \/E  and  a  parameter  vector  0  =  (e,  {bk}K}  where  eT  is  a  time  delay  expressed  as  a 
fraction  £  of  the  pulse  width  T.  The  set  {bk}K  contain  the  signal’s  symbols. 


x(t) 


x(t  ^  f  bk  Yk(t)  for  0  <  t  <  KT 

|  0  otherwise 

bk  e  {ei2nm/M  |  m  =  0, 1  1} 

for k  —  0:  (K—l) 


(20) 


Figure  4  shows  the  block  diagram  of  the  source,  channel  and  classification  process.  The 
source  generates  a  symbol  at  each  time  interval.  The  variable  K  represents  the  total  number  of 
symbols  in  the  MPSK  signal.  The  modulation  is  a  form  of  encoding  represented  by  bk  which  can 
take  M  possible  phase  values.  The  encoding  produces  the  signal  x{t)  which  is  corrupted  with 
AWGN  n{t)  and  produces  the  received  signal  r{t).  The  received  signal  is  decorrelated  using 
the  pulses  which  forms  an  orthonormal  set.  A  delay  £  is  added  at  the  input  of  the  correlator  to 
represent  the  chip  asynchronous  nature  of  the  correlator.  The  delay  takes  the  value  within  the 
range  [0, 1).  The  coefficients  {rk}K  become  a  finite  set  of  observations  which  is  processed  by 
the  classifier. 

The  CON  representation  in  (21)  allow  us  to  deal  with  a  set  of  K  random  variables  which 
are  the  decorrelated  symbols  at  the  receiver.  The  concept  is  also  applied  to  the  unnormalized 
transmit  signal  s(t)  and  the  noise  n(t)  due  to  the  linearity  of  the  process. 

yjt(e)  =  (y(t-£T):Yk(t)) 

K- 1 

*k(e)  =  (£  h  Wk{t-eT),\ffk(t)) 

k= 0  'AD 

nk{e)  =  (n(t  —  eT),  \\fk{t)) 
yk{e)  =sk(e)  +  nk(e) 
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Fig.  4  Block  Diagram  of  Preprocessing  Stage 


The  statistics  of  the  noise  coefficients  in  (22)  can  be  easily  obtained  by  calculating  (5)  from 
the  series  representation.  The  autocorrelation  formula  of  «*(£ )  is  obtained  from  the  scalar 
product  (iVo/2  S(t),  wM)- 

K- 1 

n(t)  =  £  nk(e)Vk(t) 
k= 0 


E{nfe(e)}  =  0 


(22) 


E{/i*(e)  «/(e)} 


The  output  of  the  correlator  is  normalized  so  its  variance  is  one.  The  conditional  likelihood 
is  given  by  (23).  The  SNR  y  is  assumed  be  known  prior  to  the  classification. 


p(rk{e)\Jt,bk,£) 

4(e) : 


_ _ _ e-||?*(e)||2/2 -^TfRe{rk{e)b*k}-y 

VnNo 

E 

7~Nq 


(23) 


This  model  assumes  phase  synchronization  for  simplicity.  The  conditional  likelihood  of  the 
MPSK  signal  in  (23)  is  calculated  by  assuming  independent  symbols.  The  average  likelihood 
function  is  calculated  in  the  next  step. 
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3.2.2.2  Averaging  over  Symbols 


Under  the  assumption  of  independent  symbols,  the  conditional  likelihood  function  is  given  by 
(24). 

K- 1 

X{{r}K\Je  =  M,{bk}K,e)=  Y[p{rk\M?  =  M,bk,e)  (24) 

k= 0 

The  condition  on  the  symbols  {bk}K  is  eliminated  by  taking  the  expectation  (4)  over  all 
possible  MPSK  symbol  values  and  results  in  a  product  of  sums  of  hyperbolic  cosine  functions 
given  by: 

7  K-lM/2-l 

X({r}K\jf^M,e)  =  e-^£^-K^Yl  E  cosh(y/2y Re{rk(e)  e~i2nm/M}).  (25) 

k=0  m= 0 


The  development  of  the  average  likelihood  in  [4]  resulted  in  a  familiar  solution  (26)  when 
comparing  hypothesis  M  against  M/2.  The  average  likelihood  over  { bk }  contains  the  power-law 
classifier  derived  from  an  Ad-Hoc  rule.  This  classifier  can  be  interpreted  as  follows:  if  a  signal 
belongs  to  class  M,  then  rk  raised  to  the  Mth  power  results  in  a  constant  value;  if  a  signal  belongs 
to  class  M/2,  then  the  operation  would  results  in  noise. 


3.2.2.3  CDMA  Model  in  Time  Domain 


A  time-domain  CDMA  signal  is  constructed  from  a  set  mutually  orthogonal  or  quasi-orthogonal 
B PS K  waveforms  pu(t )  for  u  G  (0, 1,  1}  given  by  (27). 


L- 1 


Puif)  =  E  Cl,u  Yl(t)  for  u  =  1,  - ;L 

1=0 


ck,u  e  {±1} 


(27) 


The  coefficients  cyM  are  chosen  to  ensure  the  orthogonality  of  the  set  of  waveforms.  Each 
waveform  will  serve  as  a  separable  coding  channel. 


(Pu(t),Pv(t)) 


L?  u  —  v 
a  0  w/v 


(28) 
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An  ideal  LxU  CDMA  signal  x(t)  is  defined  in  terms  of  the  sum  of  modulated  versions 
of  these  waveforms  with  a  set  of  modulation  parameters  {blljc}UxK  where  the  variable  K  now 
represents  the  total  number  of  chip  intervals  in  a  single  transmission.  Each  waveform  u  is 
normalized  with  a  factor  1  / \f~L  such  that  the  waveform  energy  is  1.  The  product  of  the  data 
vector  and  normalized  waveform  has  energy  Eu. 


x(t)=x(t,C,  { b}K/L ) 


k  puit-kLT ) 
0 


for  0  <  t  <  KT 
otherwise 


(29) 


The  channel  is  characterized  by  an  AWGN  process.  Similar  CON  pulses  to  (9)  are  used  at 
the  correlator.  Each  coefficient  has  two  subscripts  i  and  k  that  identify  the  location  of  the  ith 
pulse  within  the  kth  symbol  interval.  Each  coefficient  is  a  function  of  £  which  now  characterizes 
a  frame  asynchronous  version  of  the  CDMA. 


ri,k{£)  — 


1 


--{y(t —  eT),\j/i(t —  kLT)) 


VW2 

K/L-Iu-I 

Si,k(e)  =  (  £  £  bu,k  Puit  -  kLT  -  sT),  ¥i(t  -  kLT)) 

k= 0  u=0 

nt,k{£)  =  (■ n(t  —  eT),\i/i(t-kLT )) 

£  G  0, 1,  1 


(30) 


The  observation  and  noise  statistics  are  obtained  in  the  same  manner  as  those  in  (22). 


E(«^(e)}  =  0 


E  {MeKm(£)} 


2  W./  ®k,m 


(31) 


For  balanced  CDMA,  we  define  a  chip  signal-to-noise  ratio  (33)  in  terms  of  a  nominal  value 
of  the  energy  per  symbol,  noise  power  value  Nq  and  code  length  as: 


7c  = 


E 

NqL 


(32) 


The  chip-level  SNR  is  not  a  good  choice  for  a  parameter  because  of  its  dependency  on  the 
hypothesis.  A  better  choice  is  the  ratio  of  the  total  energy  density  {Ej /K)  and  the  noise  power 
density  value  Nq.  For  simplicity,  the  discussion  will  refer  to  this  parameter  as  the  density  ratio. 
The  formula  is  derived  assuming  that  the  total  energy  Ej  equals  the  energy  per  symbol  times 
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the  number  of  users  (E  ■  U )  multiplied  by  the  total  number  of  symbol  intervals  ( K/L ). 


Er  L 

Density  Ratio  =  — —  =  ys—  (33) 

A_/Vo  U 

All  this  information  is  sufficient  for  expressing  the  conditional  likelihood  of  the  Lx  U 
CDMA  signal  per  symbol  frame. 

X(rk\J?  =  {L1U},C,b,£)  = 

1  -Tf(e)rk(£)/2+y/TfrRelrf(e)}sk(£)-Ycsl(e)sk(e) 

(. kNq)l! 2 

4(e)  =  {^(e)}/=0:i-1  (34) 

^(e)  =  {^(e)}/=0:i-1 
4(o)  =  { £  cLl,bH_ky  Q;L-' 

u= 0 


This  theoretical  development  summarizes  the  existing  prior  knowledge  on  classification  of 
CDMA  signals  that  was  derived  from  MPSK  classifiers.  The  next  section  will  discuss  the  steps 
necessary  for  developing  a  compact  form  of  the  average  likelihood  (35)  after  computing  the 
expectation  over  all  the  unknowns. 


X(rk\J#’)  =  Ef{ECo 


,0^0,1  ,  —  ,CL-l,U-\ 


Uk{m\Jr  =  {L:U},C,bk,e)}}}  (35) 
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4  Modulation  Classification  Assumptions 
and  Procedures 


This  section  discusses  the  strategy  for  finding  a  simplified  form  of  the  CDMA  likelihood 
function.  The  first  approach  consists  on  solving  the  average  likelihood  for  the  simplest  case 
of  CDMA,  i.e.,  M’  =  {2,2}  and  then  extend  the  solution  to  code  lengths  of  3,  4  and  higher 
code  lengths.  For  a  hypothesis  —  {2, 2},  the  averaging  process  will  produce  2LU+L  =  64 
exponential  terms  that  would  require  simplification.  This  research  found  that  a  pattern  can  be 
obtained  with  the  assistance  of  a  symbolic  algebra  algorithm  shown  in  Appendix  A.  Finding 
a  rigorous  mathematical  proof  followed  such  discovery  and  made  it  very  simple  to  develop. 
The  success  of  this  approach  was  validated  using  ROC  curves.  Such  curves  ensured  that  the 
classification  of  simple  CDMA  signal  can  be  achieved  with  the  proposed  procedure. 


4.1  Uniformly  Distributed  Codes  and  Data 


The  first  choice  for  averaging  over  the  unknown  codes  assumesthat  each  code  coefficient  of  a 
spreading  matrix  C  is  a  uniformly  distributed  (36)  binary  antipodal  symbol.  Each  coefficient 
is  assumed  to  be  completely  independent  from  all  the  other  coefficients.  A  problem  with  this 
assumption  is  the  implication  that  any  random  code  could  be  used  for  CDMA,  which  is  a 
wrong  assumption.  Although  code  coefficients  appear  to  be  random,  a  random  selection  of  code 
coefficients  does  not  necessarily  produce  a  CDMA  code  with  its  characteristic  low  correlation 
properties. 


P(C) 

P(b ) 


1 

W 

1 

w 


(36) 


Under  the  assumptions  of  perfect  symbol  synchronization  (e  =  0),  balanced  energy  (E  = 
constant)  and  full  load  (L  =  U)  CDMA,  the  conditional  likelihood  of  one  symbol  is  given  by 
(37).  If  e  is  unknown,  then  the  probability  over  the  variable  can  be  assumed  to  be  uniform: 

m = i/i. 


l(n\M-={L,U},C,b)  = 

1  _  (37) 

j—^exp{-rk{e)Hrk{e)/ 2  +  ^/2ycRe{rk(e)H }sk(e)  -  ycfk  (e)?fc(e)) 

Several  key  propositions  and  definition  were  established  for  averaging  over  C  and  b.  First, 
we  would  deal  with  the  expectation  over  the  data  vector  and  average  over  U  unknown  users. 
Second,  we  would  deal  with  the  expectation  of  the  code  matrix  by  splitting  the  sum  over  C  in 
three  terms  that  provide  simplification. 
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Proposition  4.1  The  average  of  a  function  f(C  b)  over  C  G  (±l}ix  11 
by: 


and  b  E  {±1}U  is  given 
(38) 


Proof  4.1  Defining  q,j  =  Cjjbj  and  substituting  in  the  original  equation  gives  the  desired 
result: 


1 

w 


C  h 


=  2H7EyE/(e  i) 

2  b 

^E/(gT)^Ei 

2  fo 


=  ^E/(ci)- 


(39) 


The  next  step  requires  averaging  the  likelihood  over  the  code  matrix.  The  summation  over 
all  possible  combinations  of  is  split  in  two  summation  terms.  The  first  one  is  the  summation 
over  a  set  of  matrices  §g.  A  second  summation  is  a  set  of  amplitude  vectors  a  that  have  a 
specific  construction.  Vectors  a  have  non-negative  coefficients  a,  and  can  take  odd  or  even 
values.  The  vector  represents  all  the  possible  amplitude  levels  (40)  that  can  be  generated  from 
the  product  C  b. 

0,2,4 for  U  even 
1,3,5 ,...,£/  for  U  odd 

The  space  of  all  possible  code  matrices  will  be  partitioned  into  cells  defined  by  the  set  §3. 
Definition  4.1  The  set  §3  is  defined  in  terms  of  the  amplitude  vectors  a  as: 

=  \ C  |  C  G  (±l}LxC/  and  a,-  = 


u-\  ) 

I  cu  (41) 

j= 0  J 


In  order  to  obtain  from  a  all  the  possible  values  generated  by 
necessary  to  multiply  the  amplitude  vector  by  a  diagonal  matrix 
shown  in  (42). 

G  =  diag(g) 

g€{±  1}L 


the  product  Cb  it  would  be 
G  that  restores  the  signs  as 

(42) 
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The  average  over  C  and  b  is  replaced  by  the  average  over  the  product  of  G  and  a.  A  new 
parameter  p  can  be  interpreted  as  the  probability  of  the  sign  vector  g. 


l{rk\Jt?  =  {L,U})  = 


1 


1 


Y  nY  Y  p~7»rk/2+^2YcRe{r»  GC\}-Yc\\GC\\\2 


P  = 


1 

¥ 


(43) 


The  averaging  over  all  g  is  accomplished  by  using  Proposition  2. 

Proposition  4.2  The  average  of  the  function  exp(gTa)  over  g  e  {±1}  is  given  by: 


L—  1 


¥LeS  3  =  Tlcosh( 


an 


(44) 


i=  0 


Proof  4.2  This  proposition  is  proven  by  induction.  Consider  vectors  g \  =  [go,  ■■■■,gL-\)T , 
a\  =  [a0,...,aL-  i]r,  g2  =  [g0,...,gL-hgL]T  and  a2  =  [a0,  ...,aL-i:aL]r  with  gj  e  {±1}.  The 
proposition  should  hold  true  when  increasing  the  dimensionality  of  g\  to  g2. 

1  -t  ~  L~l 

If  ¥^eSl  “]  —  J^[ coshfaf)  is  true, 


g  i 


i= 0 


1 


then  ^ e§2  °2  ~  Y\cos¥ai)  LS  true- 

gi  i= 0 

That  is: 


2l+i 


g  2 


g  1 

1 


gL 


=  ¥Y^e§l  UcosKaL ) 


g  1 


L 

j  ~[  cosh  ( 
i= 0 


an 


(45) 


The  average  likelihood  in  (43)  is  expressed  as  a  summation  of  hyperbolic  cosine  function, 
similar  to  the  likelihood  of  MPSK  signals. 


Hh\je  =  {l,u})  = 


dfrk/2-yc\\a\\2 


L- 1 


f\  cosh( \JlycRe{r*ik}ai)  (46) 

7=0 
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Equation  (46)  was  found  to  be  identical  to  the  solution  provided  by  the  symbolic  algebra 
algorithm  presented  in  Appendix  A  by  setting  a  precision  parameter  /3  =  0.  The  implementation 
offers  reduced  performance  when  classifying  between  BPSK  and  2x2  CDMA.  The  perfor¬ 
mance  of  the  classifier  in  Figure  5)  has  a  weaker  performance  when  compared  to  a  classifier 
that  was  developed  using  a  non-uniform  probability  as  shown  in  Figure  7.  The  equation  (46) 
is  difficult  to  implement  due  to  the  large  number  of  terms  generated  by  all  the  possible  com¬ 
binations  that  a.  Each  vector  coefficient  can  take  U  / 2  values  according  to  (40)  and  there  are 
L  coefficients,  which  translate  to  ( U /2)L  different  amplitude  vectors.  For  a  4  x  4  CDMA,  the 
equation  is  extremely  complex  and  the  implementation  produced  numerical  problems.  A  more 
practical  approach  can  be  obtained  from  weighting  the  CDMA  matrices  according  to  their  low 
correlation  properties. 


4.2  Averaging  over  a  Non-Uniform  Probability  of  the  Code 


When  considering  potential  approaches  for  classifying  CDMA  matrices,  some  consideration 
must  be  given  to  the  code  design.  CDMA  matrices  are  designed  in  a  way  that  they  exhibit  low 
cross  correlation  properties  between  their  column  vectors.  The  Total  Squared  Correlation  [14] 
defined  in  (47)  is  a  parameter  that  measures  the  correlation  between  column  vectors  of  the  code 
matrix. 

L-1L-1 


'=0  7=0 


(47) 


C  —  {ci,j} 


j=0:L-l 


A  modification  of  this  metric  (48)  will  be  used  for  constructing  a  weight  referred  as  the 
probability  of  a  CDMA  code. 


L-1L-1 


- j\ 

i= 0  7=0 

-r  /^ii2  ji 


t(C)  =  \\C  C\\p  —  L“  >  0 


(48) 


Only  full  rank  matrices  will  be  considered  when  designing  a  CDMA  classifier.  If  the 
condition  of  fully-loaded  (L  =  U)  is  not  enforced,  then  it  would  be  possible  to  find  matrices 
such  that  the  TSC  is  minimum;  however,  such  matrices  would  not  qualify  as  CDMA  codes 
because  extending  the  matrices  to  a  full  rank  code  matrix  gives  no  guarantee  of  achieving  a 
minimum  TSC.  Therefore,  any  classification  of  underloaded  CDMA  signals  must  be  constrained 
to  full  rank  code  matrices  that  are  highly  uncorrelated. 

The  weighted  average  can  be  seen  as  a  filtering  process  where  some  code  matrices  are 
emphasized  over  others.  Averaging  the  likelihood  over  a  uniform  probability  of  codes  means 
that  no  particular  code  matrix  is  preferred  over  others.  This  assessment  is  incorrect  because 
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the  likelihood  in  (46)  considers  codes  that  do  not  exhibit  low  correlation  properties.  The  idea 
of  using  the  TSC  norm  for  constructing  a  probability  of  code  was  developed  with  significant 
success.  Norms  such  as  the  Maximum  Squared  Correlation  (MSC)  may  be  used;  however,  for 
the  purpose  of  achieving  correct  classification,  the  usage  of  the  TSC  proved  to  be  sufficient. 

A  desirable  property  of  the  TSC  is  that  it  is  invariant  to  permutations  of  rows  and  columns. 
This  is  easily  verified  by  expressing  the  TSC  as  a  Frobenius  Norm,  which  is  also  invariant  to 
permutations  and  arbitrary  rotations.  This  invariance  allows  us  to  remove  the  energy  term  of  the 
likelihood  out  of  the  summation  over  the  sign  vector  g. 

4.2.1  Design  of  Code  Matrices  using  TSC 

The  design  of  orthogonal  matrices  for  CDMA  is  a  vast  and  complex  field  of  study  out  of 
the  scope  of  our  research  goals.  A  brief  overview  of  the  design  theory  of  orthogonal  will 
be  provided  to  provide  the  idea  of  how  the  code  matrices  were  selected  in  our  classification 
problems.  The  discussion  starts  with  the  definition  of  Hadamard  matrices. 

Definition  4.2  Hadamard  matrix  H  of  size  Lx  L  is  a  matrix  of  elements  \hjj\  =  1  that  satisfies 
the  following  conditions: 

HT  H  —  LILxL  (49) 

Hadamard  matrices  are  constrained  to  (±l}Lxi  in  the  generation  of  CDMA  signals  using 
BPSK  symbols.  There  are  many  way  to  construct  Hadamard  matrices,  but  the  most  common 
one  is  the  so  called  Sylvester’s  construction. 


Definition  4.3  The  Sylvester’s  construction  is  given  by  the  following  recursive  formula  using 
the  Kronecker  product  [15 ]: 


H2 


1  1 
1  -1 


(50) 


Hjk  —  H2®  H2k- 1 


Hadamard  matrices  have  been  a  subject  of  study  of  mathematicians  for  years.  Several 
code  constructions  such  as  Paley  and  Williamson  constructions  exist  at  the  present  time.  An 
interesting  conjecture  states  that  Hadamard  matrices  can  be  designed  for  code  lengths  divisible 
by  4.  Unfortunately,  there  is  no  theorem  that  can  prove  the  existence  for  any  arbitrary  code 
length  divisible  by  four  [16]. 

Conjecture  4.1  A  real  Hadamard  matrix  with  must  be  of  code  length  (order)  one,  two  or  a 
multiple  of  four.  [17] 

The  following  theorem  is  a  modification  of  the  relationship  between  the  Frobenius  Norm 
and  the  Kronecker  product  [15]. 
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Theorem  4.3  The  Total  Squared  Correlation  of  a  matrix  generated  from  a  Kronecker  product 
is  given  by:  [15] 

T{Ha®Hb)  =  T(Ha)z{Hb)+L2az(Hb)+L2bz(Ha)  (51) 

Equation  (51)  leads  to  the  following  corollary. 

Corollary  4.4  If  the  TSC  of  two  matrices  are  zero,  then  the  TSC  of  the  Kronecker  product  must 
also  be  zero. 


4.2.2  Probability  of  the  Code  Matrix 


The  proposed  probability  of  a  code  matrix  is  a  function  of  the  TSC  such  that  the  weight  for 
small  values  of  TSC  are  high  while  the  weight  for  high  TSC  values  are  low.  The  first  choice  is 
to  construct  the  exponentially  decaying  probability  function  shown  in  (52). 


nc) = 

IT  =  £w(C)c~/3t(c) 


w(C) 


1  for  t(C)  ^  'Emin 
0  otherwise 


(52) 


This  probability  resembles  a  sampled  version  of  an  exponential  probability  function  with  a 
precision  parameter  /3  that  controls  the  selection  of  code  matrices.  This  definition  seems  to  be  a 
unique  contribution  of  this  effort,  because  this  concept  does  not  appear  the  after  a  review  of  the 
technical  literature.  Other  probabilities  can  be  derived  using  similar  concepts,  but  the  interest  in 
this  particular  one  is  the  preservation  of  the  exponential  nature  of  the  entire  likelihood  function. 
The  probability  can  be  tailored  for  a  limited  set  of  TSC  values  least  or  equal  than  t(C)  <  Tmax. 
We  assume  that  Tmax  — »  °°. 


4.2.3  Expectation  over  a  Non-Uniform  Probability  of  Code 


The  development  of  the  likelihood  under  the  assumption  of  non-uniform  probability  of  code 
starts  with  Proposition  4.5  and  results  in  (53)  assuming  that  a  frame  synchronous  signal  is 
provided. 


X(rk\M’={L,U}) 


1 

(n  Nq)l/2 


EEpE  E  p^e 

3  T  g  CeS^T 


rfrk/2+f2^Re{Tf  GC  l}-ye||GCl[|2 

(53) 
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The  space  of  all  possible  code  matrices  will  be  partitioned  in  non-overlapping  sets  similar 
to  the  methodology  used  for  uniformly  distributed  codes.  These  sets  contain  matrices  such  that 
the  absolute  value  of  the  sum  of  the  columns  equals  some  vector  a  and  a  code  matrix  in  the  set 
has  a  TSC  value  r  according  to  54. 


Definition  4.4  The  set  T  is  defined  as: 


1 

r 

u- 1 

|c  C e  {±l}Lxt/,a;  = 

I  CU 

7=0 

and  t(C)  =  T 


(54) 


The  summation  over  the  code  coefficients  is  split  into  three  summation  terms,  the  summation 
over  g,  C  G  §3)T,  and  r  as  shown  in  (53).  Then,  we  proceed  to  rearrange  the  summations  as 
follows: 


X{?k\M’  =  {L,U}) 


e-'fV  2 
(n  Nq)l/2 


e^pi|2e  e 

ci  T 


e-p  m 

w 


Pi 


es/WcRe{> f  Ga} 


(55) 


Averaging  over  G  is  performed  by  swapping  the  diagonal  terms  in  the  matrix  with  the  observa¬ 
tion  vector  r  in  the  correlation  term. 


rT  G  a  — gT  diagir)  a  (56) 

Proposition  4.2  can  be  applied  to  obtain  the  average  likelihood  function  in  the  form  of  products 
of  hyperbolic  cosine  functions  given  by  (57). 

L-l 

e-7cll°ll  a(d,/3)  Y\cosh(s/2yc  Re{r*k}  at)  (57) 
i= 0 

_  _  e-pm 

<*&P)  =  I  I  — nr~  (58) 

t  ce 

This  equation  can  be  divided  in  four  terms:  the  energy  of  the  received  signal,  the  energy  of 
the  model,  a  term  that  depends  on  the  code  and  a  product  of  hyperbolic  cosine  functions  that 
depends  on  the  correlation  between  the  received  signal  and  the  model. 

4.2.4  Analytical  Average  Likelihood 

The  empirical  likelihood  obtained  from  the  symbolic-algebra  algorithm  is  consistent  with  the 
derived  likelihood  (57)  and  (58)  for  code  lengths  between  1  through  4.  The  generality  of  the 
development  give  no  reason  to  suspect  that  this  average  likelihood  formula  would  not  hold  for 
higher  code  lengths. 


X(h\^  =  {L,U})  = 


(k  No)1/2 
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4.2.4.1  BPSK  versus  2x2  CDMA 


The  likelihood  (57)  of  a  single  symbol  interval  with  M'  =  {L,  U }  provides  an  exact  solution 
(60)  applicable  to  both:  uniformly  distributed  code  /3  =  0  and  the  non-uniform  probability 
of  the  code  /3  ^  0.  It  also  provides  the  formula  for  BPSK  signals  by  just  setting  L  =  1.  A 
binary  classification  between  BPSK  and  2x2  CDMA  requires  comparing  pairs  independent 
BPSK  symbols  against  single  2x2  CDMA  symbols.  We  obtain  the  formula  for  two  BPSK 
independent  symbols. 

Kin}2 =  {1, 1})  =  e~Yc  cosh(y/2yc Re{ro}) -e~Yc  cosh{y/lyc Re{r{\)  (59) 

The  exact  likelihood  of  a  2  x  2  CDMA  provides  means  to  explore  the  effect  of  the  precision 
parameter  on  the  overall  likelihood  function. 

\(rk\Jf?  =  {2,2})  =  1  m(e~sl5  +  e~4jc  cosh(2y/2yc  Re{r0jk})+ 

2(1+ e  (60) 

e~Ayc  cosh(2y/2yc  Re{rXik})  +  e~*P~Ayc  cosh(2^/2yc  Re{ro,k})cosh(2^/2yc  Re{rljk})) 

Figures  5  through  7  show  the  Receiver  Operating  Characteristic  (ROC)  Curves  for  a  LRT 
for  2  x  2  CDMA  vs.  BPSK.  These  curves  were  based  on  2000  tests,  each  test  uses  a  sequence  of 
128  chips,  and  plotted  for  different  total-SNR  (total  energy  over  No).  Each  graph  is  for  a  single 
value  of  J3.  The  best  performance  is  obtained  when  /3  — »  °°.  Due  to  the  exponential  probability 
(52)  the  coefficients  a  decays  quickly  for  low  TSC  values. 
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Probability  of  Detection 


ROC  2x2  CDMA  Vs.  BPSK,  128  chips/test,  2000  tests,  (E=0 


Fig.  5  ROC  =  {2,2}  versus  ^  =  {1,1},  j3  =0 
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ROC  2x2  CDMA  Vs.  BPSK,  128  chips/test,  2000  tests,  p  =  l 


Fig.  6  ROC  Mi  =  {2,2}  versus  Mx  =  {1, 1},  /3  =  1 
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Probability  of  Detection 


ROC  2x2  CDMA  Vs.  BPSK,  128  chips/test,  2000  tests,  p  =  °o 


Fig.  7  ROC  M’x  =  {2,2}  versus  M’x  =  {1,1},  j8  =  °o 
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4.3  Simplification  of  the  Average  Likelihood 

The  likelihood  of  a  CDMA  can  be  further  decomposed  in  two  terms:  those  that  quickly  decay 
to  zero  and  those  that  persist  when  /3  — >  The  dependency  of  each  a(a.p)  coefficient  for  a 
2x2  CDMA  is  shown  in  Figure  8  and  9.  In  this  specific  case,  a([0,0]r,/3)  and  a([2,2]r,/3) 
converge  to  zero  while  a([0, 2]r,/3)  and  a([2,0]r,/3)  converge  to  1/2. 


Log[a,([0,0]./l-»ixi)] 

a 


Fig.  8  2x2  CDMA  Vanishing  Coefficients  a  versus  /3 ,  a  G  {[0,0]r,  [2,2]r} 


The  computation  of  these  coefficients  become  extremely  difficult  even  for  simple  cases  like 
a  2  x  2  CDMA.  A  logical  approach  must  consider  the  limit  as  /3  — >  <*>  using  the  definition  of  a 
provided  in  58. 


Definition  4.5  A  feature  vector  a*  =  {a*}  for  a  given  hypothesis  MJ  —  {L.U}  is  defined  in 
terms  of  minimum  TSC  matrices  C*  such  that: 


T(C)  >  T(C*)  =  Ttmin 


U- 1 


7=0 


(61) 


Proposition  4.5  The  limit  of  oc(d.fi))  as  approaches  infinity  depends  on  the  feature  vectors 
as  follows: 


lim  a  (a  3) 

P — ^°° 


|S5 


Y*ae{a*} 


for  af  a* 
for  a  —  a* 


(62) 
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a([0,2]^-*  oo) 


a 


Fig.  9  2x2  CDMA  Non- Vanishing  Coefficients  a  versus  /3,  a  e  {[0, 2]r,  [2,0]r} 


Proof  4.3  From  a  set  of  all  possible  vectors  {aj}  and  all  possible  TSC  values  { T, }  with 
0  <  To  <  T i<  Tf+i  for  M'  —  {L,  U },  we  define  d,  j  —  |§g  T.|  and  dj  —  Y,jdij-  The  coefficients 
can  be  expressed  as  a  ratio  of  exponential  functions. 


lim  a(a,/3) 


dj,ke  Ptk 

Lidie-fo 


(63) 


Then,  we  multiply  by  e^T°  on  both  the  numerator  and  denominator.  Any  term  with  T/.  >  To 
vanishes  as  /3  — >■  °°  leaving  only  the  terms  where  T^  =  Tq. 


lim  a  (a,  J3) 

j8— »°o 


di^e 

Eid,-e-^T'~To) 


(64) 


Therefore,  the  coefficients  depend  on  the  feature  vector  set  {a*}. 


lim  a  (a,  1 3) 

P — 


0 

doj 

do 


for 

for  a  —  a* 


(65) 


By  using  an  infinite  precision,  many  terms  vanishes  from  the  average  likelihood  function 
leaving  only  the  terms  that  are  contributions  of  the  lowest  TSC  code  matrix  and  allowing  to 
perform  the  classification  in  a  computationally  efficient  way.  The  question  remains  whether 

34 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


computing  the  likelihood  coefficients  would  be  affordable  after  this  simplification.  The  Mat- 
lab  algorithm  shown  in  Appendix  B  was  implemented  for  computing  for  small  code  lengths. 
Unfortunately,  the  code  becomes  useless  for  processing  code  lengths  greater  than  4.  A  compu¬ 
tationally  efficient  approach  for  estimating  the  coefficients  can  be  constructed  from  interpreting 
the  meaning  of  ce(a,°°). 


4.3.1  Interpretation  of  the  CDMA  Likelihood  Coefficients 


The  values  of  cx(a.  fi)  can  be  calculated  exactly  for  code  lengths  of  2  as  follows.  First,  the 
feature  vectors  a  are  identified. 


=  [0.2] 

=  [2,0]r 


The  code  matrices  of  interest  in  the  sets  T  are: 


Sa^  =  (Ctt  Sa*t  =  {C2*} 


c?= 


i  i 
i  -i 


<2  = 


(67) 


With  each  set  of  cardinality  1,  it  is  easy  to  calculate  the  values  of  the  likelihood  coefficients  as 
follows: 


a(ci\,oo) 


a(<7£,°°) 


Pg*,ol  _  1 

|^a|,o|  T  |§a*,o|  2 

!§a?,ol  _  1 

l§3*,o|  +  |§a2,ol  ^ 


(68) 


A  closer  look  at  (68)  reveals  that  the  likelihood  coefficients  a  are  a  ratio  of  the  occurrence 
of  particular  CDMA  vectors  abs(C  b)  over  the  total  number  of  feature  vectors  that  exists  for 
given  hypothesis  M’  =  (L,  U}.  A  more  practical  way  of  calculating  the  coefficients  comes  from 


(69). 


a  (a* 


count (abs(C *  b)  =  a* 

7P 


(69) 


This  limit  is  an  approximation  because  for  a  given  code  length,  there  may  exist  several 
orthogonal  matrices  C*  that  are  not  related  by  a  permutation  of  columns  and  rows.  These  code 
matrices  are  referred  as  non-equivalent  matrices  [18].  Unfortunately,  despite  the  numerous  code 
constructions  that  exist  in  the  literature,  there  is  no  specific  theorem  that  predicts  the  number 
of  non-equivalent  matrices  for  a  specific  code  length.  It  is  known  that  for  code  lengths  of 
L  =  4,8,16, 32, 64  and  128  there  are  1, 1 , 5, 3, 60  and  487  non-equivalent  matrices  respectively 
[19].  This  limits  the  construction  of  the  average  likelihood  to  our  knowledge  on  codes.  It  was 
found  that  extracting  the  feature  vectors  of  each  of  the  five  16  x  16  CDMA  non-equivalent 
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matrices  resulted  in  the  same  values  of  a *  with  slightly  different  coefficients  values.  An 
assumption  that  non-equivalent  Hadamard  matrices  generate  the  same  set  of  feature  vector  will 
be  use  in  our  classification  problems.  The  algorithm  for  generating  the  likelihood  coefficients 
and  the  feature  vectors  is  provided  in  Appendix  D. 

The  following  procedure  summarizes  the  construction  of  the  simplified  likelihood  function 
given  a  hypothesis  {L,  U }.  The  procedure  is  applied  to  find  the  likelihood  function  of  a  4  x  4 
CDMA. 

1.  Find  C*,  a  code  matrix  with  minimum  TSC; 

2.  Extract  the  feature  vector  from  C*  b  by  varying  b\ 

3.  Compute  the  likelihood  coefficients  using  (69); 

4.  Using  the  coefficients  of  the  feature  vector,  construct  the  product  of  hyperbolic  cosines  in 
(57). 

The  4  x  4  Hadamard  matrix  is  given  by  the  Sylvester’s  Construction  (50).  There  are  two 
equivalent  matrices  of  interest  that  produce  two  important  feature  vectors.  All  other  feature 
vectors  can  be  derived  from  the  permutation  of  these  vectors. 

1111]  [-1111 

1-1  1  -1  1-1  1  1 
1  1  —  1  —  1  2  1  1-1  1 
1-1-1  lj  [  1  1  1-1 

The  4x4  case  has  two  feature  vectors  which  can  be  obtained  from  adding  the  columns 
of  the  previous  matrices.  The  likelihood  coefficients  do  not  change  when  a  a*  is  permuted. 
This  invariance  of  a  is  due  to  the  exchangeability  property  [20]  of  the  model  that  makes  the 
likelihood  invariant  to  row  permutations.  Because  the  transmission  model  admits  CDMA  code 
matrix  C  with  any  arbitrary  permutation  of  rows,  the  final  likelihood  does  not  depend  on  the 
order  of  the  vector  coefficients  of  r. 

The  cardinality  of  the  sets  is  computed  as  follows.  For  C\,  the  rows  that  add  up  to  zero 
can  be  permuted  in  3!  ways.  Also  there  are  24  ways  of  changing  the  sign.  The  value  |Sj*;o| 
is  3!  ■  24  =  96.  For  C4,  the  rows  that  add  up  to  2  can  be  permuted  in  4!  ways.  Also  there  are 
24  ways  of  changing  the  sign.  The  value  |§>^ j0|  is  4!  ■  24  =  384.  At  the  time  of  computing  the 
coefficients  a,  one  needs  to  account  for  4  possible  permutations  the  a*.  The  total  number  of 
Hadamard  matrices  for  a  4  x  4  code  matrix  is  the  sum  of  the  cardinalities  of  these  sets,  that  is 
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768.  This  result  is  in  agreement  with  the  experimental  computations. 


a{=[4,0,0,0]r 

a$  =  [i,i,i,i]r 


«(«2,°°) 


l§aTol 

4|S^,ol  +  |S^,o 

l§W,ol 

4|§a|,o|  +  l^a*,0 


1 

8 

1 

2 


(71) 


Finally,  the  4x4  average  likelihood  for  a  single  symbol  is  constructed  by  forming  the 
hyperbolic  cosine  products. 


e-(rkSk)/ 2  /  3i  _ 

=  {4,4})  =  ^kN^l/2  (  E  rfco5/r(4v/2^Re{r/]fe}) 

+  ^  n  e  7cco5/7  ( 2  x/2^  Re(r/l } )  ] 

2  i= 0  / 


(72) 


As  the  code  length  and  number  of  user  grow,  the  average  likelihood  becomes  a  complicated 
expression  and  further  simplification  is  needed  to  avoid  numerical  errors  or  any  excessive 
computational  burden. 

4.4  Discussion  and  Results 

Figure  10  through  1 1  show  ROC  curves  for  several  simple  cases  of  spreading  codes.  The  length 
of  the  signal  under  test  is  expressed  in  number  of  chips.  A  test  sample  of  128  and  192  chips 
was  chosen  for  code  lengths  multiple  of  2  or  3  respectively.  The  plots  are  based  on  2000  tests. 
The  total  SNR  is  the  total  energy  Et  divided  by  the  noise  power  Nq.  Few  observations  can 
be  made  from  these  simulations.  The  conclusion  is  that  the  average  likelihood  classifier  for 
CDMA  is  capable  of  classifying  small  code  lengths.  Figure  14  shows  that  it  is  possible  to 
distinguish  between  BPSK  and  CDMA  with  no  surprise  because  the  proposed  development 
follows  a  standard  procedure  that  can  be  extended  to  other  hypothesis  such  as  MPSK  or  QAM 
modulation. 

4.4.1  Average  Likelihood  for  the  Underloaded  Case 

The  underloaded  CDMA  model  assumes  that  the  data  vector  b  is  given  by  (73).  The  equation 
introduces  new  variables  A,  that  represent  whether  the  user  i  is  active  (1)  or  inactive  (0).  Defining 
a  scenario  of  U  <  L  users  can  be  expressed  in  L!/((L  —  U)\U\)  different  permutations,  but  all 
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Probability  of  Detection 


ROC  4x4  CDMA  Vs.  3x3,  192  chips/test,  2000  tests,  p=°o 


Probability  of  False  Alarm 


Fig.  10  ROC  Curve  for  =  {4,4}  vs.  {3,3}  CDMA 
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ROC  2x2  CDMA  Vs.  BPSK,  128  chips/test,  2000  tests,  (3  =  °° 


Probability  of  False  Alarm 


Fig.  1 1  ROC  Curve  for  ^  =  {2, 2}  CDMA  vs.  BPSK 
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these  results  are  equivalent  if  one  introduces  an  arbitrary  permutation  matrix  g?  between  Q 
and  b  which  leaves  the  same  model  if  Qg?  redehned  as  Q' .  It  is  not  necessary  to  average  over 
all  possible  combinations  of  A,,  because  its  average  will  be  equivalent  to  the  likelihood  of  any 
arbitrary  column  permutation. 


b  =  [b0Ao,biAi,...,bLAL]T 
A,-  e  {0, 1} 


(73) 


From  the  point  of  view  of  (57)  the  form  of  the  average  likelihood  should  remain  the  same, 
except  for  the  definitions  of  the  feature  vectors  (69).  They  must  be  calculated  assuming  that 
L  —  U  arbitrary  users  are  inactive. 

An  interesting  case  of  underloaded  CDMA  is  shown  in  Figure  12  where  a  4  x  2  CDMA 
fails  to  be  differentiated  from  a  2  x  2  CDMA.  The  behavior  of  the  classifier  is  explained  by 
understanding  the  Sylvester’s  Construction  that  was  use  to  generate  the  code.  A  4  x  2  CDMA 
can  be  seen  as  two  appended  CDMA  vectors  from  a  2  x  2  matrix.  Other  underloaded  cases  are 
shown  in  Figures  13  and  14.  Both  figures  validates  the  possibility  of  classifying  underloaded 
CDMA  for  small  code  lengths. 

Another  interesting  construction  is  the  case  of  L  x  1  CDMA.  The  Lx  1  CDMA  likelihood 
function  reduces  to  the  likelihood  of  L  independent  BPSK  symbols;  therefore,  it  will  be 
impossible  to  differentiate  between  those  two  classes  given  the  provided  CDMA  model.  In  such 
cases,  the  classification  would  have  to  rely  on  the  cyclostationary  concepts. 

The  next  section  provides  a  simplification  of  the  likelihood  function  as  well  as  an  extension 
to  much  larger  code  length  values. 
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ROC  4x2  CDMA  Vs.  2x2,  128  chips/test,  2000  tests,  p  =  °o 


Fig.  12  ROC  Curve  for  =  {4,2}  vs.  {2,2}  CDMA 
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Probability  of  Detection 


ROC  3x2  CDMA  Vs.  2x2,  192  chips/test,  2000  tests,  (3=<x> 


Fig.  13  ROC  Curve  for  =  {3,2}  vs.  {2,2}  CDMA 
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ROC  4x2  CDMA  Vs.  3x2,  192  chips/test,  2000  tests,  p  =  °o 


Probability  of  False  Alarm 


Fig.  14  ROC  Curve  for  =  {4, 2}  vs.  {3, 2}  CDMA 
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5  Method  Simplification,  Assumptions  and 
Procedures 


The  goal  of  this  section  is  to  approximate  of  the  average  likelihood  function  developed  in  the 
previous  section.  The  approximation  is  needed  for  processing  high  values  of  code  length  that 
could  quickly  produce  numerical  problems  and  computational  inefficiencies  due  to  the  large 
number  of  variables.  The  approximation  is  achieved  by  either  reducing  the  set  of  the  feature 
vectors,  using  Taylor  series  and  taking  advantage  of  the  properties  of  the  likelihood  function 
which  is  an  even  and  exchangeable  function. 

5.1  Revisiting  the  CDMA  Model 

The  CDMA  model  will  be  expressed  in  terms  of  the  matrix  equation  shown  in  (74).  For 
simplicity,  we  will  assume  that  all  the  variable  Q,  g  and  b  are  real  and  binary  antipodal.  The 
system  is  assumed  to  be  balanced,  with  energy  per  symbol  E.  The  code  matrix  is  an  L  x  U 
matrix.  The  noise  nk  is  now  a  real  Gaussian  random  vector  with  zero  mean  and  correlation 
matrix  No/2 1.  The  observation  rk  —  yk/ \J No / 2  is  also  defined  as  a  real  vector. 

}'k  = 

G  =  diag(g) 


-GQbk  +  nk 


(74) 


These  variables  are  averaged  using  the  following  weights  for  g,  b  and  Q. 

1 

P  = 


2L 


(75) 


P(Q)  = 


-jSt  (Q) 


Using  the  key  propositions  and  definitions  discussed  in  section  4,  we  find  that  our  likelihood 
model  is  given  by  equations  (76)  through  (78). 


L- 1 


h(rk\Jf  =  {L,U})  =  - —  n  £<Urf||fl|l~a(d,/3  )Y[c°sh(vS2yc  Re{ritk}  dj)  (76) 

(nNo)LG  s  f=o 


«(s,/3)  =  E  E 


r  CeSs 


-pm 

w 


(77) 
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0 


for  a  /  d* 


lim  a(a,/3) 

P — y°° 


is, 


d.Tm 


I  So,- 


for  a  =  a* 


(78) 


5.2  Classification  Using  a  Single  Feature  Vector 

A  distinctive  characteristic  of  the  Sylvester’s  construction  (79)  is  having  a  row  and  column  of 
ones.  A  feature  vector  derived  from  the  product  C  1  (79)  provides  a  simple  way  to  classify 
CDMA  signals  because  it  reduces  the  complexity  of  equation  (76). 

ak  =  [L,0,0,...0]r 

,1  J  (79) 

a(a  ,oo)  =p3 * 


The  average  likelihood  function  of  such  scheme  reduces  to  the  following  expression: 

_ ?H -p  / 2  L—\ 

A  (rk\Jf  =  {L,  U})=  *nN^L/2  e~YcL  P3*  E  co5/i(L  \f2yc  r*uk)  (80) 

The  construction  of  the  likelihood  ratio  between  hypotheses  =  {2L,  2L}  and  = 
{L,L}  is  a  simple  ratio  of  hyperbolic  cosine  functions.  The  terms  containing  the  energy  of  the 
vectors  and  No  are  cancelled  out.  The  formula  requires  to  give  a  special  attention  should  be 
given  to  the  format  of  the  observation  as  either  2L  x  2L  or  L  x  L  CDMA. 

m\2L)  _ «, 

M't-'t+lW  L^oCosh(^/2foLr'lOZlfIoCosh(s/TycoLr,ik+1)  h„ 

Figure  15  through  17  show  the  ROC  curves  for  detecting  CDMA  of  various  code  lengths 
using  the  provided  simplification.  The  plots  are  based  on  2000  tests,  with  128  chips/test,  using 
frame  asynchronous  symbols  and  plotted  for  several  Total-SNR  values.  For  small  code  lengths, 
the  feature  vector  (79)  appears  to  be  relevant  in  the  classification  of  classes,  but  as  the  code 
length  increases,  the  performance  of  the  classifier  is  degraded  under  the  same  amount  of  samples. 
This  degradation  occurs  presumably  for  two  reasons:  1.  classifying  higher  code  lengths  at  the 
same  total-SNR  requires  more  samples  and  2.  This  feature  (79)  becomes  less  significant  as  the 
code  length  increases. 

The  CDMA  classifier  represented  by  the  rule  in  (81)  can  be  interpreted  as  follows.  Within  a 
single  frame  of  L  chips,  the  classifier  must  find  a  single  spike  of  amplitude  proportional  to  L 
and  probability  p3*.  This  rule  is  only  valid  under  the  assumption  of  fully-loaded,  small  code 
lengths  CDMA. 
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Probability  of  Detection 


ROC  4x4  CDMA  Vs.  2x2,  128  chips/test,  2000  tests,  p=°o 


Probability  of  False  Alarm 


Fig.  15  Single  Feature  Detection:  ROC  {4,4}  versus  {2,2} 
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ROC  16x16  CDMA  Vs.  8x8,  128  chips/test,  2000  tests,  p  =  °o 


Probability  of  False  Alarm 


Fig.  16  Single  Feature  Detection:  ROC  {16, 16}  versus  {8,8} 


47 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


Probability  of  Detection 


ROC  32x32  CDMA  Vs.  16x16,  128  chips/test,  2000  tests,  p=°o 


Probability  of  False  Alarm 


Fig.  17  Single  Feature  Detection:  ROC  {32,32}  versus  {16, 16} 
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5.3  Taylor’s  Approximation 


The  first  consideration  in  approximating  the  average  likelihood  is  to  select  the  proper  sample  size. 
A  sample  size  of  K  chips  must  be  divisible  by  the  value  of  the  code  length  of  all  the  hypotheses, 
so  each  hypothesis  contains  an  integer  number  of  symbols  K/L.  Under  this  assumption,  the 

terms  Ilfio  1  e  '/k'rk^2  and  jTlfio  '  (kNo)l/2  are  common  factors  to  each  hypothesis  and  can 
be  neglected  if  necessary.  For  convenience,  the  energy  term  of  the  observation  is  kept. 

K/L- 1  K/L- 1  L- 1 

Y[  A(r* \jj?)=  n  £  fx{a,oo)eM~~7c^~  \\cosh(^2yc  riMcii)  (82) 

k= 0  k= 0  3e{3*}  i= 0 


The  product  of  the  hyperbolic  cosine  functions  can  be  approximated  using  (83).  The 
likelihood  expression  results  in  exponential  terms  that  depends  on  the  absolute  value  of  each 
observation  r,^.  In  both  approximations,  there  is  a  e±]  term  that  can  be  neglected  because  it  is 
a  common  factor  to  all  the  likelihood  functions. 


L— 1 

IT 

(=0 


Jog(cosh(-s/2y/  a,)) 


?Ef=o  VWc  abs[ruk)  a,— 1 


for  yc  >  1 
for  yc  <C  1 


(83) 


The  approximation  for  low  values  of  yc  did  not  have  a  significant  improvement  on  the 
classification,  so  only  the  approximation  for  yc  A>  1  was  considered  regardless  of  the  value  of 
Yc-  Combining  (82)  and  (83)  under  the  assumption  of  yc  1  result  in: 


K/L- 1 

II 


k=  0 


K/L- 1 

e-\\abs(rk)-y/2%a\\2/2-log(a(3,°°)) 

k= 0 


n  i 


(84) 


Using  the  log-likelihood  instead  of  the  likelihood  provides  further  simplification.  The  log- 
sum  of  exponents  in  (84)  can  be  replaced  by  choosing  the  maximum  value  of  all  the  exponents 
in  the  sum. 

/og(^e£0/)  ~  max(oJi) 

1  (85) 


(O,  =  ~\\abs(rk)  -  ^/2Ycdi\\2 /2  +  log(a(dhoo)) 


With  || a || 2  =  L  U  and  1 1 iy  1 1 2  constant,  the  maximum  value  of  CO,  depends  on  the  correlation 
of  a  and  absir *.)  and  the  coefficient  a(a,°°).  Not  much  simplification  can  be  applied  to  the 
likelihood  coefficient  other  than  investigating  its  behavior  through  actual  calculations.  This 
behavior  is  shown  in  Figure  18  and  shows  the  values  of  a(d,°°)  bounded  by  the  red  and  the 
blue  graphs  as  the  code  length  increases.  Our  approach  ignores  the  effects  of  these  coefficients 
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because  bounds  are  difficult  to  assess.  The  approach  also  relying  exclusively  on  the  correlation 
term  for  approximating  the  likelihood  function. 


CDMA  Likelihood  Coefficients  Boundaries 


Fig.  18  Boundaries  of  Minimum  and  Maximum  Coefficients  vs.  log(L) 

A  given  feature  vector  a*  and  all  its  permutations  qualify  as  features  for  classification.  This 
fact  presents  a  problem  because  the  number  of  permutations  becomes  large  for  higher  code 
lengths  causing  the  algorithm  to  be  computationally  inefficient  when  computing  (95).  For 
example,  ct  =  [4,4,0,0,0,0,0,0]r  is  a  feature  vector  under  hypothesis  —  {8,8},  but  also 
all  its  8 !/((8  —  2)  !2!)  =  28  permutations.  In  this  case,  there  are  two  distinct  vector  coefficients. 
In  order  to  deal  with  this  problem,  we  take  advantage  of  a  Computer  Sciences  theorem  known 
as  the  Maximum  Scalar  Product. 
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Proposition  5.1  Ifx  —  {x,}L  and  y  =  {xj}L,  and  8?  is  an  arbitrary  Lx  L  permutation  matrix, 
then  the  maximum  scalar  product  between  x  and  88  y  is  given  by: 


L—  1 

max{xT  &>y)=Y<  xo(i)  ya(i) 
i= 0 

xa(i)  —  '':ct(z'+  1 ) 

yjz(i)  —  yn{i+\) 


(86) 


where  o(i)  and  7i(i)  represent  permutations  of  the  it  l  element  with  i  e  {0, 1,2,  ...,L  —  2}. 


Proof  5.1  Let’s  assume  that  the  sorted  vectors  are  given  by  x(  a , 
claim  that: 


We  define  c[kukl_M  as: 


L- 1 

c[kuk2,...,kn\  =  52  Xi 

(=0 

i^ki 


TTien,  it  can  be  shown  that: 


(44  and  y{7t)  =  We 


(87) 


(88) 


+44  +cihM  > 442  +44  +c[8M 

is  true  since 

4t  04!  -  42 ) + ,*2]  >  4,  (42  -  4 ) + %  ,fe] 

For  distinct  indexes  k\  and  kj,  ify'kl  >  4’  t/ro  inequality  above  is  true  because: 

4,  >  42 
4,  -42 

Ify’ki  >  4 ,  //ir  inequality  also  remains  true  because: 

44-1 +4242 +c[fci,Jt2]  -  4i42  +44  +c[ll2] 

and 

4  >  4, 

4  -4,  >o. 


(89) 

(90) 

(91) 

(92) 

(93) 
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Once  a  permutation  of  two  arbitrary  indexes  is  verified,  we  can  proceed  with  a  permutation  of 
two  arbitrary  indexes  £3  and  k4  contained  in  tfi^kly 


Ai  y'kx  +x'k2y'k 2 +xkiy'k3  +xkfy'u +c[kUk2MM] 


>  Xkfy'k2  +4-3^4  +^4^3  +Xk2y,k1  +  C[k1,k2,k3,k4] 

(94) 


The  inequality  also  holds  whether  xki  >  xk4  or  yk}  >  yk4  or  xk4  >  xki  oryk4  >  yky  The  process 
of  permuting  indexes  is  continued  until  the  original  sorted  sequences  become  two  arbitrary 
permuted  ones:  x  =  {xt}L  and  y  =  {yi}L-  ■ 


The  use  of  Proposition  5.1  allows  estimating  the  maximum  value  of  (95).  The  process  is 
depicted  in  Figure  19.  A  sequence  of  decorrelated  samples  from  the  CON  decorrelator  retains 
only  the  absolute  value  of  the  real  part.  The  samples  are  formatted  as  an  L  x  1  vector.  The 
vector  components  are  sorted  and  sent  into  a  classifier  that  computes  the  likelihood  (95)  for  a 
given  hypothesis  (L,  U }.  The  same  process  is  repeated  for  multiple  hypotheses.  Finally,  a  MAP 
classifier  determine  the  class  by  computing  the  maximum  log-likelihood  of  all  the  hypotheses. 

max(£\ffk\a*)  =  RT  A 

i 

R  =  sort(\rk\)  ^ 

A  —  sort(af) 


Fig.  19  Block  Diagram  of  the  Computation  of  the  CDMA  Log-Likelihood 


5.3.1  Extraction  of  Features  and  Coefficients 

The  feature  vectors  for  a  given  hypothesis  must  be  extracted  and  stored  in  a  database  prior  to 
the  classification.  The  extraction  process  is  obtained  by  using  the  algorithm  in  Appendix  D. 
Sorting  the  vectors  and  computing  the  rate  of  occurrence  reduces  the  complexity  that  comes 
from  dealing  with  different  permutations  of  the  feature  vector. 
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In  the  case  of  perfectly  orthogonal  codes,  the  feature  vectors  could  be  computed  by  solving 
the  Diophantine  equations  given  by  (96).  The  features  can  take  values  of  a,  e  {2n}n=0'2'L  and 
fn  is  the  number  of  occurrences  of  a;  as  given  by  (40),  in  other  words,  \\a\\2  —  LU.  The  solution 
can  be  obtained  by  using  the  Lenstra-Lenstra-Lovasz  lattice  basis  reduction  algorithm  [21]. 

u-i 

E  (2 n)2fn=LU 

i= 0 

(96) 

i= 0 

0  <  fn  <  L,  fn  G  Z* 

5.4  Averaging  over  Unbalanced  Energy 

The  expectation  over  the  energies  of  different  users  can  be  computed  under  the  following 
assumptions.  A  nominal  energy  E  per  symbol  is  assigned  to  each  user.  The  energy  variation 
from  the  nominal  value  is  represented  as  a  continuous  random  variable  A,  such  that  the  energy 
per  symbol  of  the  user  i  is  given  by: 


Ei=E(  1+A,;)2. 


(97) 


The  model  of  the  CDMA  is  constructed  such  that  the  likelihood  can  be  separated  into  two 
terms:  a  contribution  of  a  balanced  CDMA  and  a  deviation  8. 

y  —  \J E  7 LCb  +  <5  +n  (98) 


Under  the  assumption  that  A  is  a  diagonal  matrix  with  uniform  or  normally  distributed 
coefficients,  the  vector  8  =  y/E/L  C  D  b  can  be  represented  as  Gaussian  random  vector  (using 
the  Central  Limit  Theorem)  with  zero  mean  and  covariance  matrix  a2  /.  The  expected  value  of 
the  conditional  likelihood  on  8  is  integrated  over  the  interval  (— °°,  +<*>)  as  shown  in  (99).  For 
simplicity,  the  integral  consider  only  the  terms  that  depends  on  8. 


E5.{<?  7c8‘  cosh(y/2ycri%kai)} 


-tc  <y2 

e  \/ 1+  2  yc  it2  - 

=cosh(y/ 2  Ycnjcdi ) 
Vl+2  YcG2 


(99) 


Finally,  the  average  likelihood  function  of  the  unbalance  CDMA  is  the  product  of  the 
average  likelihood  of  the  balance  CDMA  and  a  factor  %  given  by  (100).  The  addition  of  any 
additional  users  with  unbalanced  energy  decreases  the  likelihood  and  degrades  the  accuracy  of 


53 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


the  classifier. 


-Uyc  o2 
£  y/l  +  2  yc  a2 


(100) 


XiYcO2) 


(1  +  2  yc  c2)u/2 


5.5  Discussion  and  Results 


Our  simulations  are  based  on  the  code  presented  in  Appendices  C  to  F.  Each  plot  assumes 
that  the  code  can  appear  in  any  form  of  permutation,  the  observation  is  frame  asynchronous 
with  £  e  [0,L  —  1].  The  tests  consist  of  binary  classifications  using  a  MAP  approach  using 
log-likelihood  functions.  The  probability  of  the  classes  is  uniform,  which  implies  that  the  same 
number  of  test  samples  is  the  same  for  each  class.  A  samples  of  K  chips  is  used  to  compute  the 
likelihood  function.  The  sample  contains  K/L  symbols. 

Figure  20  and  21  show  the  accuracy  versus  total-SNR.  It  was  found  that  a  energy  density  to 
noise  density  ratio,  referred  as  density  ratio,  serves  as  a  better  parameter  for  all  code  lengths. 
Usage  of  the  density  ratio  parameter  used  in  Figures  22  through  26  reveals  the  existence  of  a 
minimum  required  density  ratio  value  for  accurate  classification  between  classes. 

The  plots  are  based  on  2000  tests  with  chips  lengths  that  varies  between  26  to  2 16  depending 
on  the  hypothesis.  The  chip  length  must  ensure  that  there  is  enough  number  of  symbols  in  a 
class  for  make  a  classification.  The  required  chip  length  appears  to  be  a  function  of  the  symbol 
length  of  the  highest  class. 

The  simplified  likelihood  function  does  not  assume  that  the  codes  are  perfectly  orthogonal. 
As  long  as  the  code  matrix  achieves  a  minimum  TSC,  the  equation  should  remain  valid.  The 
likelihood  function  for  CDMA  can  be  used  to  classify  CDMA  from  other  signals,  given  that  the 
log-likelihood  functions  of  the  other  signals  are  provided. 

The  simulations  shown  in  Figures  27  to  31  demonstrate  that  the  concept  can  be  extended 
to  underloaded  CDMA  cases.  We  may  find  examples  where  the  classification  is  ambiguous  in 
underloaded  CDMA  cases.  A  number  of  active  users  less  that  the  code  length,  U  <L  can  be 
confused  with  a  L/2  x  U  CDMA  especially  when  the  signals  are  constructed  using  (50). 

A  particular  case  of  interest  is  the  classification  between  orthogonal  and  non-orthogonal 
CDMA  matrices.  Figure  32  shows  the  classification  between  classes  {63,63}  and  {64,64}. 
The  feature  vectors  were  obtained  from  a  64  x  64  Sylvester  Hadamard  matrix  on  both  cases. 
The  feature  vectors  of  the  63  x  63  CDMA  were  obtained  by  removing  an  arbitrary  column  and 
row  of  the  orthogonal  matrix.  At  the  moment  of  testing  classes,  the  Hadamard  with  63  rows 
was  substituted  with  a  Gold  Code  matrix  made  from  the  following  polynomials:  x8  +  x  +  1  and 
x8  +x7  +x2  +X1  +  1.  Because  the  Gold  Code  shares  similar  features,  it  was  possible  to  classify 
between  both  classes. 
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Classification:  Full  Loaded,  Balanced,  {1,1}  vs.  {2,2} 
CDMA,  2000  tests,  128  chips/test 


Fig.  20  Classification  {1,1}  versus  {2,2}  Using  Symbol  SNR 

5.5.1  Execution  Times 

The  simulations  were  conducted  in  a  cloud  using  Matlab  parallel  programming.  A  single  test 
was  assigned  to  one  of  48  available  nodes.  The  estimated  time  of  execution  is  per  test  is  given 
in  Table  4.  The  execution  times  depend  on  several  variables,  such  as  the  data  transfer  times 
between  nodes,  the  number  of  chips/test  and  the  total  number  of  feature  vectors. 
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Classification:  Full  Loaded,  Balanced,  {2,2}  vs.  {4,4} 
CDMA,  2000  tests 
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Fig.  21  Classification  {2,2}  versus  {4,4}  Using  Symbol  SNR 


Table  4  Execution  Times  per  Test 


Hypothesis  ^ 

Hypothesis 

chips/test 

features 

Execution  Time  (sec.) 

{1,1} 

{2.2} 

128 

2 

3.3 

{2,2} 

{4,4} 

128 

3 

3.6 

{4,4} 

{8,8} 

128 

5 

4.0 

{8,8} 

{16,16} 

512 

21 

30.3 

{16,16} 

{32,32} 

2048 

56 

39.5 

{63,63} 

{64,64} 

4032 

98337 

15268 
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Classification  Full  Loaded,  Balanced,  {1,1}  vs.  {2,2} 
CDMA,  2000  tests,  K=128  chips/test 


Fig.  22  Classification:  {1,1}  versus  {2,2}  Using  Density  Ratio 
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Accuracy 


Classification:  Full  Loaded,  Balanced,  {2,2}  vs.  {4,4} 
CDMA,  2000  tests 


Density  Ratio  =  Et/No/K  (dB) 
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Fig.  23  Classification:  {2,2}  versus  {4,4}  Using  Density  Ratio 
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Classification:  Full  Loaded,  Balanced,  {4,4}  vs.  {8,8} 


Fig.  24  Classification  {4,4}  versus  {8, 8} 
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Accuracy 


Classification:  Full  Loaded,  Balanced,  {8,8}  vs.  {16,16} 
CDMA,  2000  tests 


Fig.  25  Classification:  {8,8}  versus  {16, 16} 
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Classification:  Full  Loaded,  Balanced,  {16,16}  vs.  {32,32} 
CDMA,  2000  tests 


Fig.  26  Classification:  {16, 16}  versus  {32,32} 
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Accuracy 


Classification:  Underloaded,  Loaded,  Balanced  CDMA 
2000  tests,  K=192  chips/test 


Fig.  27  Classification  of  Underloaded  CDMA 
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Accuracy 


Classification:  Underloaded,  Loaded,  Balanced  CDMA 
2000  tests,  K=192  chips/test 


Fig.  28  Classification  of  Underloaded  CDMA  for  Code  Length  L 
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Accuracy 


Classification:  Underloaded,  Loaded,  Balanced  CDMA 
2000  tests,  K=64  chips/test 
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Fig.  29  Classification  of  Underloaded  CDMA  for  Code  Length  L 


64 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


Accuracy 


Classification:  Underloaded,  Loaded,  Balanced  CDMA 
2000  tests,  K=128  chips/test 


Fig.  30  Classification  of  Underloaded  CDMA  for  Code  Length  L  — 
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Accuracy 


Classification:  Underloaded,  Loaded,  Balanced  CDMA 
2000  tests,  K=256  chips/test 
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Fig.  3 1  Classification  of  Underloaded  CDMA  for  Code  Length  L  —  32 


66 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


Classification:  Full  Loaded,  Balanced  63-Gold  Code  vs. 
64  Hadamard  CDMA,  2000  tests,  K=2~12  chips/test 


Fig.  32  Classification  of  {63,63}-Gold  Code  versus  {64,64}-Hadamard 
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Table  5  CDMA  Study  Cases 


CDMA  Schemes 

BPSK-data 

QPSK-data 

BPSK-code 

(spreading) 

Type  1 

Ce{± \}LxU 
be{±l}Uxl 

Type  3 

C6{± \}LxU 
be{±l,±i}Uxl 

QPSK-code 

(spreading) 

Type  2 

Ce{±l±i}LxU 

be{±l}Uxl 

Type  4 

Ce  {±l±i}LxU 
be{±l1±i}Uxl 

6  Extension  of  the  Average  Likelihood  Method 
Assumptions  and  Procedures 


The  development  of  the  likelihood  function  for  LxU  CDMA  was  presented  in  the  previous 
section  and  now  it  is  time  to  extend  the  algorithm  to  ah  types  of  CDMA  [3]  generated  from  a 
combinations  of  BPSK  and  QPSK  symbols  as  shown  in  Table  5.  An  alternative  definition  can 
be  obtained  by  applying  90  degree  rotations  to  types  2-4,  but  because  C  and  b  are  unknowns, 
the  selection  of  the  model  is  a  preference  rather  than  a  necessity.  A  90  degree  rotation  only 
means  that  the  real  and  imaginary  parts  of  the  model  can  be  exchanged  without  affecting  the 
modulation  type. 

The  new  developments  start  from  reformulating  the  models  in  terms  of  real,  binary  antipodal 
matrices.  The  three  terms  of  the  likelihood  function  will  be  computed:  the  correlation  term  rT  s, 
the  energy  term  sT  s  of  the  conditional  likelihood  and  the  covariance  matrix  of  each  CDMA 
model  nT  n.  The  observation  and  noise  vectors  are  given  by  a  combination  of  the  real  and 
imaginary  parts  of  the  CDMA  vectors. 


6.1  General  CDMA  Model  for  Complex  Types 


The  conversion  of  a  complex  model  to  a  real  model  is  possible  when  we  use  the  homomorphic 
transformation  (101)  of  complex  matrices  in  CLxL  to  M2Lx2L  real  matrices. 


rR  +  ir,  ->• 


7r 

n 


nR  +  ini  — y 


nR 

ni 


Qr  +  iQi 


Qr 

Qi 


-Qi 

Qr 


(101) 
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For  complex  codes,  we  define  a  new  version  of  the  TSC  in  (102).  The  variable  Lq  represents 
the  effective  code  length,  i.e.,  the  actual  code  length  of  the  complex  model  times  a  factor  of 
2.  The  probability  of  the  code  matrix  (48)  in  terms  of  the  TSC  retains  its  form.  Also,  we  may 
consider  extending  the  sets  §3  and  §3  T  to  include  complex  codes  in  the  form  of  block  matrices. 

T(G)  =  WQT  ■  QWf  ~ Lq  ■  rank(Q)  (102) 

Definition  6.1  The  subset  Sg+^_  of  matrices  Q  —  Q  \  +  iQ  -  with  Q+  and  Q  G  {±1  jLxlJ  is 
defined  as: 

^a+,3-  —  {2  I  —  2+1  and  a _  =  Q  1  with  ( a+)i  >  0  and  («_),-  >  0  }.  (103) 

Definition  6.2  The  subset  §3  3_jT  ls  defined  as: 

$a+j-,T  =  (2  !  2  e  S3+53_  W  t(Q)  =  t}.  (104) 


The  next  two  properties  will  allow  simplifying  the  CDMA  model  in  terms  of  binary  antipodal 
matrices.  Both  properties  are  used  in  the  development  the  average  likelihood  for  types  2  and  4. 


X 


/  / 
-7  / 


jr1 


x  ■ 


Qr 
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-2/ 

Qr 


■A’1 


Qr  —Qi 
Qi  Qr 


6.2  Type  3  CDMA  Average  Likelihood 


(105) 

(106) 


The  type  3  CDMA  Average  Likelihood  is  very  simple.  We  start  with  the  block  matrix  model 
in  (107).  There  are  two  orthogonal  equations  of  interest,  one  for  yn  and  the  other  for  yj.  By 
treating  the  real  and  imaginary  part  as  independent  equations,  the  problem  can  be  divided  in 
two  sets  of  orthogonal  Type  1  CDMA  under  the  assumption  of  uncorrelated  noise.  Matrices 
Gr\  and  Gr2  are  two  diagonal,  binary-antipodal  matrices  that  has  the  possible  sign  variations 
of  the  real  and  imaginary  parts  of  the  transmit  CDMA  signal. 


yn 

yi 


Gri  0 

o 

O) 

bR 
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to 
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—  Gri  Qr  ijR  +  h.R 
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(107) 
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The  average  likelihood  function  for  type  3  can  be  approximated  by  the  product  of  the 
likelihood  in  (107)-I  and  (107)-II.  Although  the  approximation  assumes  that  Gam  /  Gr 2,  the 
simulations  show  that  the  assumption  is  reasonable  for  detecting  type  3  without  increasing  the 
complexity  of  the  average  likelihood  function. 

^TlijR  +  iri\J4?)  —  ^ti(/r\J^)  ■  kTi{ri\Jrt?)  (108) 

An  empirical  calculation  of  the  likelihood  using  symbolic  algebra  reveals  that  the  true 
likelihood  is  a  product  of  hyperbolic  cosine  terms  that  depend  on  r^aR,  rfcij,  rJaR  and  rjaj. 
The  exact  equation  ensures  that  G’ai  =  Gr2-  This  relationship  is  not  significant  at  the  time  of 
classifying  CDMA  in  terms  of  the  code  length  and  number  of  users  for  type  3,  so  it  can  be 
relaxed  to  Gr\  7^  Gri. 

6.2.1  Discussion  and  Results 

The  experiments  assume  that  the  code  matrices  appear  in  any  form  of  column  and  row  permu¬ 
tation  and  the  correlator  process  frame  asynchronous  samples.  The  complex  noise  power  is 
Nq/2,  so  the  real  and  imaginary  components  contribute  with  half  of  the  noise  No/ A.  Type  3 
CDMA  uses  the  same  features  for  Type  1  CDMA.  The  average  likelihood  is  the  product  of  two 
orthogonal  Type  1  CDMA,  so  the  performance  is  expected  to  be  similar.  Figure  33  shows  a 
comparison  between  the  classification  of  the  classes  {2,2}  and  {4,4}.  The  accuracy  of  type  3 
has  a  similar  profile,  but  shifted  backward  by  approximately  3  dB. 

The  plots  in  Figure  35  were  based  on  a  variable  chips  lengths  that  varies  between  26  to 
2 16  depending  on  the  hypothesis.  The  chip  length  must  ensure  that  there  is  enough  number  of 
symbols  in  a  class  for  make  a  classification. 

A  characteristic  feature  of  these  plots  is  the  sudden  break  in  the  performance  of  the  classifier. 
The  minimum  density  ratio  appears  to  be  a  function  of  the  highest  code  length  hypothesis  in  the 
classifier.  Figure  34  shows  an  approximate  relationship  between  the  75  percent  accuracy  point 
and  versus  the  highest  code  length  in  a  binary  classification.  (The  data  was  interpolated  from 
the  markers  in  35  with  the  code  length  is  shown  in  the  legend.)  An  interpolation  on  the  graph 
suggests  that  a  {16, 16}  versus  {24,24}  classifier  requires  a  minimum  density  ratio  of  6.14  dB. 
The  experimental  value  agrees  with  this  prediction. 
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Classification  Type  1  and  3,  K=128  chips/test,  2000  tests 
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Fig.  33  Classification  of  {1, 1}  versus  {2, 2}  for  Types  1  and  3 
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Minimum  Density  Ratio 


Break  Point  Density  Ratio  (-75%  accuracy)  vs.  code  length 


Fig.  34  Classifier  Breaking  Points  for  Type  3 
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Fig.  35  Classification  of  {L,L} 


versus 


jL/2,L/2] 


for  log2(L)  between  2  and  7 
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6.3  Type  2  CDMA  Average  Likelihood 


The  average  likelihood  of  a  type  2  CDMA  can  also  be  expressed  in  terms  of  the  type  1 
likelihood  under  simplifying  assumptions.  The  model  in  (109)  resembles  a  type  1  CDMA  with 
the  hypothesis  rtf’  —  (2L,  U }  where  U  represents  either  the  number  of  variables  in  the  real  or 
imaginary  part.  Its  block  matrix  G  consists  of  four  diagonal  matrices  with  symbols  ( Gr and 

(Gi)i,i  e  {0,±1}. 
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The  complex  code  matrices  in  this  study  were  derived  from  Hadamard  codes.  The  construc¬ 
tion  was  made  by  selecting  a  Hadamard  matrix  and  multiplying  it  by  a  permutation  involutory 
matrix  to  form  the  imaginary  part.  A  90  degree  rotation  in  the  complex  plane  might  be  required 
in  order  to  express  the  code  according  to  the  definition  in  Table  (5). 

Definition  6.3  A  permutation  involutory  matrix  has  the  following  properties: 

PP  =  I, 

P  1  =  1,  (110) 

PiJ  e  {0.1}- 

Proposition  6.1  If  H  is  a  real  Hadamard  matrix  and  P  is  a  permutation  involutory  matrix  with 
P  f  /,  then  the  code  Q  =  H  +  iP  H  is  a  complex  Hadamard  matrix. 

Proof  6.1  This  can  be  verified  by  calculating  QH  Q  —  (2 L)I.  ■ 


Based  on  the  previous  proposition,  one  can  derive  matrices  Qr  and  Qj  satisfying  the 
condition  imposed  in  Table  5.  Setting  the  sign  of  the  amplitude  vectors  to  G/i  =  Gp_  —  0, 
Gri  ^  Gr2  =  0  results  in  two  orthogonal  sets  of  equations  given  by:  (1 1 1). 


I  :  yR  = 

II  :yi  = 


Gri  Qr  bR  +  nR 
Gr2  Qi  t>R  +  n/. 


(Ill) 


Without  entering  in  details,  the  average  likelihood  approximation  of  type  2  is  identical  to  the 
type  3  already  discussed. 


Arefa  +  iri\J^)  =  InijR  +  irfiM’) 
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A  computation  using  symbolic  algebra  reveals  that  the  true  likelihood  function  of  type  2  is 
a  hyperbolic  product  of  the  terms:  r^aR,  rjaj,  rjaR  and  rjaj.  The  exact  equation  characterizes 
a  true  type  2  signals  that  satisfies  the  condition  Gr\  =  Gr2-  This  likelihood  would  be  useful  for 
differentiating  between  CDMA  types;  however,  the  relationship  is  not  essential  for  classifying 
CDMA  in  terms  of  the  code  length  and  number  of  users. 

6.3.1  Discussion  and  Results 

Figures  36  to  38  show  the  classification  of  type  2  signals  using  200  tests,  and  chip/test  that 
varies  from  28  to  216.  The  tests  assume  that  the  code  matrices  comes  in  any  form  of  column  and 
row  permutation  and  the  correlator  process  frame  asynchronous  samples.  The  complex  noise 
power  is  Nq/2,  so  the  real  and  imaginary  components  contribute  with  half  of  the  noise  No/ 4. 
The  graphs  are  similar  to  the  type  3  CDMA  classification  that  was  already  discussed. 

Figure  37  shows  several  classifiers  and  the  characteristic  breaking  points  for  fully-loaded 
balanced  CMDA.  Figures  39  to  38  shows  the  classification  of  unbalanced  CDMA  with  variances 
in  the  amplitudes  o2  =  {0,  .1,  .25}.  The  effect  on  amplitude  unbalance  on  the  accuracy  is  more 
evident  in  the  classification  of  higher  code  lengths  since  its  effect  grows  exponentially  as  a 
function  of  the  number  of  users  according  to  (100). 
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Fig.  36  Classification  of  {L,L} 


versus 


jL/2.L/2] 


for  log2(L)  between  2  and  7 
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Classification  Type  2  unbalanced  CDMA  K=1024  chips/test,  200  tests 


Energy  Density  =  Et  /  No  /  K  in  dB 


Fig.  37  Classification  of  {1, 1}  versus  {2, 2}  for  Types  1  and  3 
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Accuracy 


Classification  Type  2  unbalanced  CDMA  K=2048  chips/test,  200  tests 


Energy  Density  =  Et  /  No  /  K  in  dB 

Fig.  38  Classification  of  {1, 1}  versus  {2, 2}  for  Types  1  and  3 
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Classification  Type  2  unbalanced  CDMA  K=1024  chips/test,  200  tests 


Energy  Density  =  Et  /  No  /  K  in  dB 


Fig.  39  Classification  of  {1, 1}  versus  {2, 2}  for  Types  1  and  3 
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6.4  Type  4  CDMA  Average  Likelihood 


The  type  4  CDMA  model  can  be  viewed  as  a  type  1  CDMA  with  twice  the  code  length  and 
number  of  users.  The  model  is  subject  to  code  matrices  that  follows  the  structure  given  by 
(113). 
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If  the  assumptions  Gri  —  Gr2,  Gn  =  Gn  are  enforced,  then  the  sets  SgSjg/iT  need  to  be 
considered  in  the  development  of  the  likelihood.  The  feature  vectors  Ur  =  QRbR  —  Qpi  and 
aj  =  QibR  +  QrI>i  have  positive  coefficients  by  definition.  The  product  over  the  block  matrix 
G  would  restore  the  signature  ±1  to  each  feature  vector,  which  is  a  0  or  a  180  degree  rotation 
in  the  complex  plane.  We  should  also  allow  these  block  matrices  to  implement  a  90  and  -90 
degree  rotations  because  such  rotations  do  not  change  the  modulation  type.  In  all  these  cases, 
the  diagonal  matrices  Gri  and  G/i  can  take  values  (0, +  1,  — 1}  according  to  the  rotation  in 
the  complex  plane  such  that  any  combination  given  by  ±(aR)  j  ±  i(a;)j  or  ±(aj)j  ±  i(aR)j  is 
included  in  our  averaging  process. 

The  conversion  to  a  model  made  of  binary  antipodal  variables  (1 14)  is  obtained  by  multiply¬ 
ing  the  matrix  X  in  equation  (105)  and  applying  the  property  discussed  in  (106).  Vectors  \bR-,bj], 
[nR;Uj\  and  [rA> ; r/]  are  redefined  using  the  following  convention:  =  —xR  +xj,  x+  =  xr+xj 

to  produce  the  binary  antipodal  vectors  in  (115).  The  simplest  model  is  obtained  by  setting 
Gn  =  Gj2  =  0. 
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Qr  - Qi 
Qi  Qr 


(115) 


The  average  likelihood  for  type  4  CDMA  (116)  is  a  special  case  of  type  1  CDMA  with 
twice  the  code  length  and  number  of  users: 

)iTA(jR  +  ihW  =  {L,U},N0/2)  «  Xti(R\JE  =  {2L,2U},N0),  (116) 

where  R  =  [r+;r_],  N  =  [n+;H_],  E{NNT}  =  Nq.  The  real  component  of  the  noise  is  assumed 
to  satisfy  E {nRn\}  =  E {UjnJ}  —  Nq/ 2. 
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6.4.1  Alternative  Model 


We  define  g+  and  to  be  binary  antipodal  vectors  {±1  }L  satisfying  the  following  relationship: 


g+=gR  +  gi 

8 -  =  ~8r  +  8i 
Gm  =  diag(gR ) 
G/i  =diag(gi). 


(117) 


We  can  calculate  the  correlator  term  of  the  conditional  likelihood  using  (113)  with  Gr\  =  GR2 
and  G/i  =  Gn  and  take  advantage  of  the  Proposition  4.2. 
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The  following  form  of  the  correlation  provides  what  is  necessary  to  compute  the  average  over 
the  block  matrix  G. 


g+ 
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Qr  - Qi 

?+ 
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Qi  Qr 
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(119) 


Since  the  block  matrix  G  is  orthonormal,  the  energy  term  remains  the  same: 
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Qi  Qr 
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Qi  Qr 

b- 

di 

(120) 


The  difference  between  the  type  4  model  in  (1 19)  and  the  type  1  correlation  term  in  (56)  is 
the  substitution  of  diag(r)  with  the  following  matrix: 


diag{jR)  — >  D  — 


diag(-rj) 
diag(rR ) 


diagirR ) 
diag{rj ) 


(121) 


Using  a  development  a  similar  treatment  used  in  type  1  CDMA,  we  can  arrive  to  the  average 
likelihood  function  shown  in  equations  (122)  and  (123). 

-RtR/2  2L-1  _ 

?iT4{R\d^  =  {L,U})  =  J>-yc||a|1  a(a,P)  cosh(^/2yc  (122) 

(kNq)lG  Y  f=o  j 
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with  R  —  [r+;r_],  a  =  [or; a/]  and 


a(a,/3) 


=  E  E 

T  GeS^a,. 


e-^T(G) 

w 


£ 

2LNq' 


(123) 


The  same  formula  can  be  use  for  Type  2  CDMA,  which  can  be  seen  as  a  special  form  of 
type  4  CDMA  with  b[  =  0. 


6.4.2  Discussion  and  Results 


The  first  task  in  generating  simulations  was  the  extraction  of  the  feature  vectors  defined  in 
(104).  It  was  found  that  the  composite  feature  vector  A  =  [or;  a/]  has  the  same  features  for  a 
corresponding  Type  1  of  size  (2L,  2 U}.  This  finding  supports  our  approximation  of  a  L  x  U 
Type  4  CDMA  likelihood  to  a  Type  1  2L  x  2U  CDMA.  The  experiments  also  assume  that  the 
code  matrices  comes  in  any  form  of  column  and  row  permutation  and  the  correlator  process 
frame  asynchronous  samples.  The  noise  power  of  the  model  is  Nq/2  and  consists  of  the 
contribution  of  the  real  and  imaginary  components  under  the  assumption  of  uncorrelated  noise. 
Only  the  a  short  range  of  density  ratios  is  considered  due  to  the  long  execution  time. 

Figure  40  resembles  the  same  characteristic  features  of  the  previously  discussed  classifiers. 
It  seems  that  all  complex  types  of  CDMA  can  be  viewed  as  a  different  way  of  encoding  binary 
antipodal  symbols.  The  average  likelihood  under  the  unbalanced  energy  assumption  degrades 
the  performance  of  the  classifier  as  expected. 
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Classification  Type  4  CDMA,  K=8192  chips/test,  200  tests 


Fig.  40  Classification  of  {L,L}  versus  {L/2,L/2}  for  log2{L )  between  2  and  7 
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7  Conclusion 


The  knowledge  on  modulation  classification  applied  to  CDMA  was  limited  the  use  of  cyclosta¬ 
tionary  features  in  neural  network  prior  to  this  research.  Prior  state  of  the  art  was  limited  to 
single  user  CDMA  signals  and  relied  entirely  on  the  cyclostationary  properties  to  detect  code 
lengths  by  using  a  feature  based  neural  network  with  multiple  layers.  Such  theory  cannot  be 
easily  extended  to  the  classification  of  multiuser  CDMA  signals. 

A  generalized  likelihood  approach  was  considered  as  a  potential  method  for  classifying 
CDMA.  The  main  weakness  of  this  approach  is  the  costly  estimation  of  the  code  matrices  and 
the  data  vectors,  especially  if  one  consider  the  estimation  of  the  parameters  under  the  wrong 
hypothesis. 

Th  main  contribution  of  the  research  presented  in  this  document  was  able  to  provide  a 
classification  algorithm  for  multiuser  CDMA  using  the  average  likelihood.  Due  to  the  model 
limitations,  the  algorithm  cannot  be  applied  to  single  user  CDMA.  However,  by  developing 
decision  theoretic  approach,  we  ensure  that  the  classification  of  multiuser  CDMA  is  optimal 
by  virtue  of  the  theory.  Although  the  approximation  of  the  average  likelihood  may  be  turn  the 
classifier  in  a  suboptimal  method,  it  is  the  only  available  method  at  this  moment. 

The  construction  of  the  average  likelihood  function  was  based  on  the  following  assumptions: 

■  chip-synchronous, 

■  frame-asynchronous, 

■  arbitrary  codes  achieving  lower  bound  TSC, 

■  unknown  data  vector, 

■  unknown  permutations, 

■  fully-loaded, 

■  underloaded, 

■  balanced,  and 

■  unbalanced  energies. 

The  development  of  the  average  likelihood  for  CDMA  was  possible  when  key  concepts  such 
as  the  probability  of  the  code  matrix  and  key  properties  such  as  the  average  of  the  exponential 
functions  were  formulated.  Averaging  over  a  weight  that  depends  on  the  TSC  behaves  as  a 
filtering  process  where  matrices  with  low  TSC  are  preferred  over  all  other  matrices.  By  defining 
a  precision  parameter  /3,  we  are  able  to  control  the  width  of  the  probability  of  code  matrices. 
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Taking  the  limit  of  /3  to  infinity  allows  finding  a  simplified  likelihood  function  in  terms  of 
the  feature  vectors,  which  are  unsigned  versions  of  the  product  between  the  code  matrix  and 
the  data  vector.  This  simplification  eliminates  many  negligible  terms  that  have  no  significant 
contribution  to  the  average  likelihood  function.  Furthermore,  a  Taylor  simplification  reduces 
the  likelihood  to  a  sum  of  exponentials.  Also  log-likelihood  expression  of  a  sum  of  exponentials 
can  be  simplified  to  the  selection  of  the  maximum  value  of  dot  product  of  the  sorted  versions  of 
the  feature  vectors  a*  and  the  amplitude  of  the  observation  vectors  abs{T\). 

The  proposed  method  produced  desirable  results  for  classifying  CDMA  with  code  lengths 
between  2  and  128.  The  proposed  approximation  avoids  multiple  problems  such  as  large 
number  of  unknowns,  numerical  overflow,  and  intractability  issues  derived  from  an  increasingly 
large  number  of  row  permutations  in  the  code  matrix.  Extending  the  classifier  to  higher  code 
lengths  is  always  possible;  however,  the  implementation  will  require  longer  processing  times 
for  acquiring  the  feature  vectors  and  computing  the  average  likelihood. 

The  existence  of  a  minimum  energy-density  to  noise-density  ratio  is  a  characteristic  feature 
of  the  newly  developed  CDMA  classifiers.  The  classifier’s  performance  below  this  limit 
approaches  to  50  percent  accuracy,  with  a  quick  transition  to  100  percent  accuracy.  The  slope 
of  the  accuracy  versus  density  ratio  varies  depending  on  the  number  of  samples  per  test.  The 
curve  is  steeper  with  higher  samples  per  test.  In  machine  learning  terms,  this  could  be  mean 
that  the  slope  of  a  sigmoid  function  depends  on  the  number  of  test  samples.  In  information 
theoretic  terms,  this  behavior  could  imply  that  the  capacity  of  the  channel  has  been  exceeded  to 
a  point  where  it  is  impossible  to  identify  the  modulation  type  of  the  signal.  These  hypotheses 
would  require  verification. 

These  concepts  can  be  extended  complex  types  of  CDMA.  The  entire  complex  model  can 
be  mapped  into  a  pure  real  CDMA  system  that  has  binary  antipodal  symbols.  Types  2  and  3 
CDMA  turn  into  two  orthogonal  sets  of  type  1  CDMA.  Type  4  CDMA,  under  the  assumption  of 
perfectly  orthogonal  signals  can  be  represented  as  a  type  1  CDMA  with  twice  the  code  length 
and  twice  the  number  of  users.  Because  of  this,  the  same  average  likelihood  function  can  be 
applied  to  all  types  of  CDMA  with  minor  modifications. 

CDMA  signals  may  resemble  Gaussian  noise.  However,  framing  the  signal  at  the  proper 
code  length,  removing  the  polarity  and  sorting  the  vector  coefficients  of  the  observation  remove 
the  random  nature  of  CDMA  and  make  the  signal  completely  detectable.  In  the  case  of  one 
user,  the  model  becomes  equivalent  to  an  ordinary  BPSK  or  QPSK  signal.  In  CDMA  cases 
where  the  number  of  users  less  or  equal  than  half  of  the  code  length,  the  classification  may  fail 
when  compared  to  a  signal  with  half  of  the  code  length  due  to  way  that  Hadamard  codes  of 
code  length  L  can  be  constructed  from  codes  with  code  lengths  of  L/2. 
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7.1  Future  Work 


There  are  many  scenarios  for  detecting  CDMA  that  can  be  explored.  Several  scenarios  such  as 
multipath  or  chip  asynchronous  CDMA  can  be  studied  and  tested  in  the  field.  A  more  ambitious 
development  includes  the  estimation  of  the  code  matrix  using  a  average  likelihood  and  the 
automation  of  a  demodulation  process  for  CDMA. 
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Appendix  A  Mathematica  Code  for  Com¬ 
puting  the  Average  Likelihood 
Function 

(*  CDMA  LIKELIHOOD  FUNCTION  *) 

(*  1.  Code  Length  and  Number  of  Users  *) 

L  =  '1\U  —L\ 

(*  2.  Code  Matrix  Defined  *) 

H  =  Array  [Subscript^-,  #1  —  1 ,  #2  —  1]  &,  {L,L}]; 

(*  3.  Data  Vector  Defined  *) 

B  —  Array  [Subscript^,#—  1,K]&,  {L,  1}]; 

(*  4.  Observation  Vector  Defined  *) 

R  =  Array  [Subscript^,  #1  —  1  ,K]&,  {L,  1}]; 

(*  5.  Correlation  Term  *) 

X  -  R.H.B; 

X  =  Expand[X[[l,  1]]] 

(*  6.  TX  Signal  and  Energy  *) 

Y  =  (H.B ); 

ET  =  (y.y); 

ET  =  Expand  [ET]; 

ET  =  ET/.v_A2->l; 

ET  =  ET[[1, 1]] 

(*  7.  Total  Squared  Correlation  *) 

H2  =  #.//; 

V2  =  0; 

For[il  =  l,il  <  L,il-H-, 

For[i2  =  l,i2  <  L,i2++, 

If[il^i2, 

V2  =  V2  +  H2[[il,i2]]A2 


]; 

wl  -  Expand[V2]/.{_A2  ->•  1}/.{_A4  -►  1}; 
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FML  =  wl 


(*  8.  Normalization  Constant  of  Probability  of  the  CDMA  Code  *) 
Z  =  Join  [Flatten  [//]]; 

Iflz.^1, 

Wl  =  Exp[— j3  *  wl]; 

For[i  =  l,i  <  Length[Z],i++, 

If  [  !NumberQ[Z[[i]]] , 

rl  —  {Z[[i]\  — >  z}; 

Wl  =  Wl/.rl; 

Wl=Sum[Wl,{z,{— 1,1}}] 

]; 

], 

Wl  =Exp[— j3]; 


(*  9.  Conditional  likelihood  *) 

CL  —  Exp[Sqrt[2y]  *X  —  y*  ET  —  J3  *  FML] 

(*  1 0  .AveragingtheCond  itionalLikel  ihood *) 
AL  -  CL; 

Z  =  Join[Sqrt[l]  *  Flatten [//],  Flatten [f?]] 
For[i  =  l,i  <  Length[Z],i++, 
If[!NumberQ[Z[[i]]], 
rl  —  {Z[[i]\  — »  z}; 

AL  =  AL/.rl; 

If  [NumberQ  [Z[[i]]/.  {c_._  — >•  0}] , 

AL  =  Sum[AL,{z,{— 1, 1}}], 

AL=  l/2*Sum[AL,{z,{— 1,1}}] 

]; 

]; 

]; 


(*  11.  Expansion  of  the  Exponents  *) 

ALS  =  AL/.Exp[aaa_] >  Exp[Expand[aaa]]; 

rulel  =  {Exp  [Plus[exl _ ]  +ex2 _ *  rs_.*-]  :— >•  Exp[Plus[exl]]  *  Subscript^, s] ATimes[ex2] } ; 

ALS  =  ALS//.rulel; 

rule2  =  {Exp  [ex2 _ *rs_,^] :— >•  Subscript^,  .s]ATimes[ex2]}; 

ALS  =  ALS//.rule2; 

rule3  =  {Subscript^,  s_]ATimes[ex2 _ ] » Cosh[Times[ex2]  *  Subscript^, s,  AT]] 
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+Sinh  [Times  [ex2]  *  Subscript^,.?,  A]]}; 
ALS  =  ALS//.rule3; 

ALS  =  Expand  [ALS] 


(*  12.GENERAL  FORM  OF  THE  2x2  CDMA  likelihood  Function  *) 

ALSS  =  Expand[ALS/Wl] 

4e~*i 6  4e-4rCosh[2v/2V7r0^]  4e-4rCosh  [2y/2y/yrliK]  4«-V-8rCosh[2v^v^^]Cosh[2VVrr1^] 

$+8e~sP  +  8+Se-W  +  8+8e~8^  +  8+&e~sP 


(*  13.  Precision  to  Infinity  *) 
LF1  =  ALSS/.{j3  -4  oo} 


4rCosh 


2V2^frQ.K]  +  ie-^Cosh  \ly/2y/ynjc 
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Appendix  B  Computation  of  the  CDMA 

Average  Likelihood  Coefficients 


o. _ 

o 

%  Calculation  of  the  numerator  of  the  alpha  parameters  using  formula. 
%  A.  Vega,  Last  Modified:  October  8,  2013 
%  inputs:  a  =  vector  vector,  U  =  number  of  users 
%  output:  matrix  =  [  TSC,  occurrences  ] 

%  formatted  output  =  occurrences  *  eA  (beta*TSC) 

o, _ 

o 

function  output  =  alphaParam (a, U) 

%  Ensure  positive  coefficients 
a  =  abs (a) ; 

%  Code  length  of  the  spreading  matrix 
L  =  length ( s ) ; 

%  Calculate  possible  column  permutations  such  that  sum (cols)  >=  0 
colPerm  =  []; 

colPerm ( 1 ) . x ( 1 , : )  =  ones(l,U); 
for  kl  =  l:floor(U/2) 
x  =  ones ( 1 , U) ; 
k2  =  kl; 
while  (  k2  >  0  ) 
x(l,k2)  =  -1; 

k2  =  k2  -  1; 

end 

colPerm (kl+1 ). x  =  unique (perms (x) ,' rows ' ) 
end 

%  Display  possible  row  vectors 
if  (  0  ) 

for  k3  =  1  :  size  (a ( : ) , 1 ) 
squeeze (a (k3 ) .x(:, :) ) 
end 
end 

%  Construct  matrices  using  row  vectors 
group  =  abs (U-a) /2+1; 
elems  =  zeros (1,U); 
for  kl  =  1 : L 

elems (l,kl)  =  size (colPerm (group (kl )). x, 1 ) ; 
end 

%  This  algorithm  is  inefficient  for  L>6 
tic 

stack  =  [];  mtx  =  zeros (L,U); 

start  =  ones(l,L); 

while (  all ( start<=elems )  ) 

%  Generate  Matrix 
for  kl  =  1 : L 

mtx  (kl,  : )  =  colPerm (group (kl) )  ,x(start(l,kl) ,  :) ; 
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end 

%  Calculate  Total  Square  Correlation 

tau  =  0; 

h2  =  mtx' *mtx; 
for  kl  =  1 : L 
for  k2  =  1:L 
if (  kl  ~=  k2  ) 
tau  =  tau  +  h2(kl,k2)A2; 
end 
end 
end 

%  Store  value 
stack  =  [stack, tau]; 

%end 

%  This  algorithm  explores  all  possible  combinations  of  a  row 
%  The  sum  of  the  row  >=  0 
k2  =  1; 

while (  k2  <=  L  ) 

c  =  start(l,k2)  +  1; 
if (  c  >elems(l,k2)  ) 
if (  k2  <  L  ) 
start  ( 1 , k2 )  =  1; 

k2  =  k2  +  1; 

else 

start  ( 1 , k2 )  =  c; 

k2  =  k2  +  1 ; 

end 

else 

start ( 1 , k2 )  =  c; 

k2  =  L+l ; 

end 

end 

end 

toe 

%  Count  all  unique  values 

vals  =  unique (stack f ,' rows ')' ; 

svals  =  zeros (size (vals) ) ; 

for  kl  =  1 : length (vals ) 

svals (l,kl)  =  sum (stack==vals (kl ) ) ; 

end 

%  Display  results 

[vals 1 , svals ' *2AL/2A (sum (s==0) ) ] 
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Appendix  C  Generation  of  Type  1  CDMA 

Signal 


o. _ 

o 

%  Matlab  Code  Snippet 
%  Generation  of  CDMA  signals  Type  1 
%  Inputs  : 

%  L  =  code  length,  power  of  2,  L  >  1 
%  U  =  active  users 

%  chips  =  signal  length  in  chips,  multiple  of  L 
%  delta_var  =  variance  of  the  amplitude 
%  Outputs: 

%  s  =  CDMA  signal  samples 

o. _ 

o 

%  Generate  Code  Matrix  Type  1 
h  =  hadamard (L) /sqrt (L) ; 

%  Implement  Row  and  Column  Permutation 
r  =  randperm(L); 
c  =  randperm(L); 
h  =  h  (r,  c)  ; 
h  =  h  (  :  ,  1  :  U )  ; 

%  Generate  Data  Vector 
symbols  =  randsrc (U, chips/L) ; 

%  GenerateAmplitudevariation 
if (  deltavar  ~=  0  ) 

AA  =  diag (delta_var/U*randn (U, 1 ) ) ; 

else 

AA  =  zeros (U) ; 

end 

%  CDMA  signal 
s  =  h* (eye (U) +AA) ^symbols; 

%  reformat  as  a  stream 
s  =  reshape (s,  [size (s,  1) *size (s, 2) , 1 ] ) ; 
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Appendix  D  Feature  Extraction  Algorithm 


o_ _ 

o 

%  Matlab  Code  Snippet 
%  Feature  Vector  Extraction 
%  Inputs: 

%  L  =  code  length 
%  N  =  size  of  the  histogram 
%  Outputs: 

%  aa  =  feature  vector  array 
%  alpha  =  alpha  coefficients 

o, _ 

o 

function [  aa,  alpha  ]  =  surveyf eat (L, U, N) 

hL  =  hadamard(L); 

fv  =  zeros (L, 0 ) ; 

alpha=  zeros (1,0); 

for  k  =  1  :  N 

%  Compute  CDMA 
if (  U  ==  L  ) 

s  =  hL*randsrc (L, 1 ) ; 

else 

u  =  randperm(L); 
u  =  u  ( 1 : U)  ; 

s  =  hL ( : , u) *randsrc (U, 1 ) ; 

end 

%  Sort  entries  by  absolute  value 
ss  =  sort (abs  (s) ) ; 

%  Enter  first  element 
if ( isempty ( fv)  ) 
fv  =  [  ss  ]  ; 
alpha=  [  1  ] ; 

else 

%  search  closest  match 

done  =  0;  rindex  =  1;  cindex  =  1; 

while (  done  ==  0  ) 

if (fv (rindex, cindex)  ==  ss (rindex, 1)  ) 

%  Smart  sorting  to  save  time 
if  ( rindex  ==  L  ) 

%  found  match,  go  to  next  column 
F(l, cindex)  =  F(l, cindex)  +  1; 
done  =  1; 

elseif ( rindex<  L  ) 

%  go  next  row  and  compare 
rindex  =  rindex  +  1; 

else 

%  reached  end  of  row 
if(cindex<  size(fv,2)  ) 

%  go  to  next  column 
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cindex  =  cindex  +  1; 
rindex  =  1; 

else 

%  reached  last  entry,  append  to  end 
fv  =  [  fv,  ss  ] ; 
alpha=  [  F,  1  ] ; 
done  =  1 ; 

end 

end 

elseif ( fv ( rindex, cindex)  >ss (rindex, 1 )  ) 

%  found  bigger  entry 
if (cindex  ==  1  ) 

%  append  at  beginning 
fv  =  [  ss,  fv  ] ; 
alpha=  [  1,  alpha] ; 
done  =  1 ; 

else 

%  append  in  before  the  bigger  entry 
fv  =  [  fv ( : , 1 : (cindex-1 ) ) ,  ss,  fv (:, cindex : end)  ]; 
alpha=  [  F  ( 1 , 1 :  (cindex-1 )) ,  1,  F ( 1 , cindex : end)  ]; 

done  =  1 ; 

end 

else 

%  found  smaller  vector 
if (cindex  ==  size(fv,2)  ) 

%  no  more  vectors  found,  append  at  end 
fv  =  [  fv,  ss  ] ; 
alpha=  [  F,  1  ] ; 
done  =  1 ; 

else 

%  more  vectors  available  ,  increase  indexes 
cindex  =  cindex  +  1; 
rindex  =  1; 


end 

end 

end 

end 

end 

aa  =  fv; 

alpha  =  alpha/sum (alpha) ; 
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Appendix  E  Calculation  of  the  Average  Like¬ 
lihood 


o, _ 

o 

%  Computation  of  the  log-likelihood  per  signal 
%  Inputs: 

%  rx  =  unformatted  observation  vector  with  proper  length 
%  gam  =  total  energy  to  noise  ratio 
%  L  =  code  length 

%  U  =  active  users 

%  delta_var  =  variance  of  the  amplitude 
%  Outputs : 

%  llkh  =  total  log-likelihood 

o, _ 

o 

%  Load  feature  vectors  and  alpha  coefficients 
[aa, alpha]  =  loadfeatures (L, U) ; 

%  Number  of  possible  delays,  frame-asynchronous 
delayavg  =  L; 

%  Number  of  samples  assuming  multiple  of  L 
chips  =  length (rx); 

%  factor  for  amplitude  variations 
gammac  =  gam/chips/U; 

%  log-likelihood  per  delay  per  sample 
dist  =  zeros (chips/L, delayavg) ; 

%  factor  for  amplitude  variations 

logchi  =  unbalanced_f actor (delta_var, L, U, gammac) ; 

%  Calculate  likelihood  for  different  delays 
for  delay  =  0:delayavg-l 

%  implement  a  given  delay 
r_delay  =  circshift (rx, delay) ; 

%  For  each  sample 
cnt  =  1; 

for  kl  =  1 : L : lim 

%  Format  as  a  vector  assuming  code  length  =  L 
r_vec  =  r_delay (kl : kl+L-1 , 1 ) ; 

dist2  =  loglikelihood (r_vec, aa, gammac)  +log  (alpha)  tlogchi; 

%  feature  vector,  log  sum  of  exponents 
dist  (cnt, k2  +  l)  =  log_sum_exp (  dist2  ); 
cnt  =  cnt  +  1 ; 

end 

end 

%  Compute  likelihood  per  delay 
llkh_delay  =  sum(dist); 

%  Average  over  delays 

llkh_symbols  =  log_sum_exp (llkh_delay) -log (delayavg) ; 

%  Sum  of  log  likelihood  of  all  independent  symbols 
llkh  =  sum (lkh_symbols)  ; 
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o, _ 

o 

%  Matlab  Code  Snippet 

%  Computation  of  the  distance  term  of  the  likelihood  per  feature  vector 
%  Inputs  : 

%  r  =  unformatted  observation  vector  with  length  L 
%  gammac  =  chip  level  SNR 
%  aa  =  collection  of  feature  vectors 
%  Outputs: 

%  llkh  =  log-likelihood  per  sample 

o. _ 

o 

functiondist  =  loglikelihood (  r,  aa,  gammac  ) 

%  assuming  that  feature  vector  is  already  sorted 

%  aa  =  sort (r) ; 

rs  =  sort (abs (real (r) )) ; 

dist  =  zeros (1, size (aa, 2) ) ; 

for  k=l : size (aa, 2 ) 

dist(l,k)  =  -l/2*sum((r  -  sqrt (2*gammac) *aa ( : , k) ) . A2 ) ; 

end 
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Appendix  F  Error  Probability  Matrix  Code 


o_ _ 

o 

%  Matlab  Code  Snippet 

%  Computation  of  the  error  probability  matrix 
%  Inputs: 

%  HT  =  array  of  hypothesis,  column-wise  pair  {  L,  U  } 

%  Et  =  total  energy.  No  =  noise  power 
%  chips  =  chip  length  of  the  signal 
%  L  =  true  code  length,  U  =  true  active  users 
%  delta_var  =  variance  of  the  amplitude 
%  ntests  =  number  of  tests 
%  Outputs: 

%  confmatrix  =  error  probability  matrix 

o_ _ 

o 

%  Matlab  parallel  for  loop  goes  here 
parfor  tests  =  l:ntest 

%  Generate  a  stream  of  CDMA  signals 

s  =  randBaseBandSignal (' typel ', chips , L, U, delta_var) ; 

%  Transmitted  signal  (baseband) 

Es  =  Et/ (chips/L) /U; 
tx  =  sqrt (Es) *s; 

%  choose  a  random  delay  between  0  and  L-l  chips 
delay  =  randi ( [ 0 , L-l ] , 1 , 1 ) ; 
tx  =  circshift (tx, delay ) ; 

%  AWGN  channel 

n  =  noise (NoReal ,  length (tx)  ); 
rx  =  tx  +  n; 

%  For  all  the  hypothesis  test  in  array  HT 
for  k2  =  1: length (HT) 

Ltest  =  HT (1, k2) ; 

Utest  =  HT (2, k2) ; 

gammac  =  Et/NoReal; 

rr  =  real (rx) /sqrt (NoReal/2 ) ; 

%  Store  the  true  class  and  the  loglikelihood  in  a  structure 
lkhTest (tests) .class  =  kl; 

lkhTest (tests) . Ikh(k2)  =  cdmaLogLikelihood (rr, gamma, Ltest, Utest, delta_var) ; 

end 

end 

%  Compute  Error  Probability  Matrix  using  MAP,  in  this  case  the  same  as  LRT 
for  test  =  l:ntest 

trueClass  =  lkhTest (test ). c; 

[~, decision]  =  max ( lkhTest (test ). lkh) ; 

confmatrix (trueClass , decision)  =  confmatrix (trueClass , decision) +1 ; 

end 


99 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


Nomenclature 


Acronyms 

ALRT  .  Average  Likelihood  Ratio  Test 

AGWN  .  Additive  White  Gaussian  Noise 

BPSK  .  Binary  Phase  Shift  Keying 

CDMA  .  Code  Division  Multiple  Access 

CON  .  Complete  Orthonormal  Set 

dB  .  Decibels 

EM  .  Expectation  Maximization 

FFT  .  Fast  Fourier  Transform 

GFRT  .  Generalized  Fikelihood  Ratio  Test 

HFRT  .  Hybrid  Fikelihood  Ratio  Test 

FRT  .  Fikelihood  Ratio  Test 

MMSE  .  Minimum  Mean  Squared  Error 

MPI  .  Multiple  Programming  Interface 

MSC  .  Maximum  Squared  Correlation 

PN  .  Pseudo  noise 

QAM  .  Quadrature  Amplitude  Modulation 

QFRT  .  Quasi-Average  Fikelihood  Ratio  Test 

QPSK  .  Quaternary  Phase  Shift  Keying 

ROC  .  Receiver’s  Operating  Characteristic 

SNR  .  Signal-to-Noise  Ratio 

TSC  .  Total  Squared  Correlation 
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Symbols 


(}L 

jjLxt/ 

rl 


z* 

Re{-} 

Im{-} 

(v) 


S(t) 

(■)* 

(Y 

abs(-) 

! 

7r 

n 

C h)i 

0 

1 

0 

I 

H 

IlxL 

Qr 

Qi 

ci,j 

a 

8 

dicig(g) 

P 

E 

Ej 


set  of  L  elements  or  column  vector 
matrix  of  L  rows  and  U  columns 
real  domain  of  dimension  L 
complex  domain 

domain  of  real  non-negative  integers 

real  part  operator 

imaginary  part  operator 

Kronecker  product 

inner  product 

Frobenius  norm 

vector  norm 

Dirac  delta  function 

Kronecker  delta 

Hermitian  operator 

transpose  operator 

absolute  value 

factorial  symbol 

real  part  of  a  complex  column  vector 

imaginary  part  of  a  complex  column  vector 

ith  component  of  a  real  column  vector 

null  column  vector 

column  vector  of  1  entries 

null  matrix 

identity  matrix 

Hadamard  matrix 

Lx  L  identity  matrix 

real  matrix 

imaginary  part  of  a  matrix 
matrix  element,  row  i  and  column  j 
CDMA  amplitude  vector 
CDMA  sign  vector 

diagonal  matrix  construction  based  on  g 
Precision  parameter 
Energy  per  symbol 
Total  energy 
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Nq/2  .  Noise  power  density 

y  .  symbol-level  SNR 

yc  .  chip-level  SNR 

b  .  data  vector 

C,  Q  .  code  matrices 

s(t)  .  transmit  signal 

y(t),y  .  received  signal 

r(t),  r  .  observation 

n(t),n  .  noise  signal 

E{-}  .  expectation  operator 

&  .  Autocorrelation  matrix 

&(At)  .  Autocorrelation  function 

.  Normal  distribution 

.  hypothesis 

R  .  Observation 

L  .  code  length 

U  .  active  users 

K  .  number  of  chips 

X(r .  average  likelihood 

A  ( r |  Jf?,  C,b )  .  conditional  likelihood 
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